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I.  INTRODUCTION 


The  term  radome  was  coined  during  World  War  II  as  a 
composite  of  two  words,  radar  and  dome.  The  term  originally 
described  the  class  of  dome-shaped  structures  designed  to 
house  and  protect  radar  ancennas  on  airborne  platforms. 
Current  usage  of  the  term  has  evolved  to  encompass  any 
structure  that  protects  or  houses  an  antenna  [Ref.  1] . 

Radomes  must  satisfy  a  diverse  range  of  engineering 
requirements.  These  specifications  can  be  broadly  divided  into 
two  categories:  electrical  and  aero-mechanical.  Due  to  aero- 
mechanical  considerations  (which  are  not  addressed  in  this 
thesis)  an  ogive-shaped  radome  is  commonly  found  on  airborne 
platforms.  Chapter  II  discusses  the  geometry  for  an  ogive¬ 
shaped  body  of  revolution.  Aircraft  radomes,  especially  those 
used  at  supersonic  speeds  are  subject  to  mechanical  stress  and 
aerodynamic  heating  so  severe  that  the  electrical 
characteristics  of  the  radome  may  be  a  secondary 
consideration. 

The  goal  of  this  thesis  is  to  develop  a  computer  program 
that  provides  an  accurate  model  for  the  transmission 
characteristics  of  a  radome,  given  the  geometry  and  electrical 
characteristics  of  its  material  composition.  This  program  is 
called  LDBORMM. F,  and  is  described  in  detail  in  subsequent 
chapters. 
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A.  ELECTRICAL  CHARACTERISTICS  OF  RADOMES 


Electromagnetic  transmission  characteristics  are  a  primary 
electrical  design  problem  for  radomes.  In  this  thesis,  the 
radome  is  modeled  as  a  thin  isotropic,  homogeneous  dielectric 
shell.  The  reflection  coefficient  (T)  of  the  radome  is  given 
by 


r 


_ 3 _ 

2  Rg  +  l\ 


(1.1) 


where  Rs  is  the  surface  resistivity  of  the  radome  and  r]  is  the 
impedance  of  free  space  (120tt  ohms)  .  The  surface  resistance, 
which  can  be  complex  in  general,  is  given  by 


R 


e 


£t 

Js 


(1.2) 


where  Et  is  the  tangential  electric  field  at  the  surface  of 
the  radome  and  Js  is  the  surface  current  [Ref  2:  p.  433]. 

From  Equation  (1.1),  it  follows  that  the  reflectivity  ( r ) 
of  the  radome  increases  as  the  surface  resistance  (Rs) 
decreases.  In  the  limit,  where  R,=0  (perfect  conductor) ,  r=-l 
as  expected.  When  an  antenna  radio’  _s  through  a  radome,  a 
small  portion  of  the  main  br  ‘eflected  off  of  the  walls 
and  forms  a  lobe.  This  '  _mmonly  called  the  reflection 
lobe  or  image  lobe.  e  lobe  is  usually  broader  than  the 
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The  primary  concern  regarding  the  image  lobe  is  that  large 
discrete  reflectors  in  the  detection  field  can  become  visible 
to  the  radar  through  this  side  lobe,  and  be  interpreted  as 
moving  targets  by  the  radar  [Ref.  3:  pp.  301-311].  The  reason 


✓ 

Image  lobe  / 


In  side  lobe  levels.  V 


Figure  1.1.  Radome  Reflection  Lobe  Geometry. 

for  false  detection  of  motion  from  stationary  reflectors  is 
the  apparent  doppler  of  the  return  of  the  image  lobe  since  it 
is  at  a  different  angle  and  thus  a  different  frequency  than 
the  main  beam  lobe. 


The  ideal  radome  is  perfectly  transparent  to  the 
electromagnetic  energy  that  must  pass  through  it.  The 


The  ideal  radome  is  perfectly  transparent  to  the 
electromagnetic  energy  that  must  pass  through  it.  The 
materials  available  and  suited  for  radome  construction  in  high 
velocity  airborne  systems  are  not  transparent.  Impedance 
mismatches  (radome/air  interface)  and  radome  geometry  alter 
the  transmission  characteristics  of  the  system.  The  presence 
of  a  radome  can  affect  the  gain,  beamwidth,  sidelobe  level  and 
the  direction  of  the  boresight,  as  well  as  change  the  VSWR  and 
the  antenna  noise  temperature  [Ref.  4:  p.  265]. 

B.  METHODS  OF  ANALYSIS 

Most  analytic  models  for  radomes  are  based  on  microwave 
optics  (ray  tracing  techniques) .  These  models  assume  a  radome 
curvature  large  enough  so  that  the  surface  can  be  considered 
locally  flat.  In  addition,  multiple  reflections  and  surface 
waves  are  usually  ignored.  These  approximations  are  not  valid 
for  systems  where  the  radome  is  in  the  near  field  of  the 
antenna  or  the  curvature  of  the  radome  causes  rapid  variation 
in  the  direction  of  the  tangential  vector  along  the  surface  of 
the  radome.  Consequently,  ray  tracing  techniques  are 
inaccurate  for  the  analysis  of  high  performance  antenna  and 
radome  configurations.  Figure  1.2  illustrates  the  inaccuracies 
inherent  to  ray  tracing  techniques.  Other  analytical  methods 
are  available,  but  have  various  shortcomings  and  limitations. 
Table  1.1  gives  a  summary  of  various  analytic  methods  for  the 
solution  of  electromagnetic  scattering  problems. 
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The  electric  field  integral  equation  (EFIE)  and  the 
method  of  moments  is  employed  in  this  thesis.  The  objective  is 
to  obtain  the  solution  for  the  unknown  current  density  which 


TABLE  1.1.  SUMMARY  OF  ANALYTIC  METHODS. 


METHOD 

SHAPE  *** 

Special 

Function 

Series 

Solution 

spheres, 

infinite 

cylinders, 

spheroids 

rigorous 

formulation; separation  of 
variable 

techniques ; special 
functions 

Extended 

boundary 

condition 

method 

homogeneous 

objects 

rigorous  integral 
equation 

formulation;  harmonic 
expansion  of  Green's 
functions  bodies 

Physical 

optics 

large  planar 
surfaces 

currents  deduced  from 
incident  fields 

Geometrical 

optics 

large  surfaces 

discontinuous  field 
patterns  obtained  at 
shadow  and  reflection 
boundaries 

Geometrical , 
uniform, 
physical 
theory  of 
diffraction 

large  convex 
surfaces  with 
edges  and 
discontinuities 

generalization  of 
physical  and  geometrical 
optics 

Frock  theory 

large  smooth, 
convex  surfaces 

formulation  based  on 
principle  of  locality 

Hybrid 

methods 

large  class  of 
scatterers 

integral  equation 
formulation;  combines 
method  of  moments  and 
optic  derived  methods 

Method  of 

moments, 

finite 

element  and 

finite 

difference 

large  class  of 
bodies, 
inhomogeneous 
obj  ects 

rigorous  integral 
equation  formulation 

ray  tracing  techniques  assume 
source  rays  are  parallel  and  the 
scattering  surface  is  locally  flat. 


tangential  vectors  on  the 
surface  vary  rapidly  in 
direction. 


Figure  1.2.  Ray  Tracing  Approximations. 
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occurs  in  the  integral  equation.  The  integral  equation  is  then 
solved  using  a  numerical  procedure  called  the  method  of 
moments  (MM)  [Ref.  5:  p.  671].  When  properly  implemented,  a 
method  of  moments  solution  is  considered  rigorous;  all  of  the 
scattering  mechanisms  are  included.  This  thesis  adapts  a  code 
developed  by  J.R.  Mautz  and  R.F.  Harrington  for  scattering  by 
bodies  of  revolution. 

The  proximity  of  the  antenna  to  the  radome  requires  a  near 
field  solution  of  the  integro-dif ferential  electric  field 
equations.  Chapter  III  develops  the  analytical  solution  for 
a  radome  in  the  near  field  of  a  circular  antenna  and  the 
numerical  method  for  its  calculation. 

Chapter  IV  provides  details  of  the  program  LDBORMM . F .  This 
program  solves  for  the  far  field  radiation  pattern  of  a 
circular  aperture  transmitting  through  an  arbitrary-shaped 
dielectric  body  of  revolution.  The  EFIE  for  electromagnetic 
scattering  is  solved  numerically  for  the  system  depicted  in 
Figure  1.3. 

Chapter  V  discusses  the  testing  and  evaluation  of  the 
main  program.  The  effects  of  radome  geometry  and  dielectric 
properties  on  the  far  field  radiation  pattern  are  examined. 

The  appendices  provide  copies  of  the  source  codes  and  a 
description  of  variables  in  the  arguments  of  the  subroutines 
and  functions  contained  in  the  main  program. 
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linearly  polarized  ,  main  beam 
antenna  electronically 


scanned 


Figure  1.3  Physical  Arrangement  for  LDBORMM.F. 
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II.  OGIVE  GEOMETRY 


A.  BACKGROUND 

Forward-looking  radar  systems  in  missiles  and  aircraft  are 
by  necessity  located  in  the  nose  of  the  platform.  Aerodynamic 
considerations  dictate  an  ogive-shaped  radome  for  most 
applications.  The  basic  problem  is  illustrated  in  Figure  1.3. 
In  this  thesis,  the  subroutine  OGIVE  provides  the  physical 
dimensions  of  the  surface,  where  the  integro-dif ferential 
field  equation  solutions  are  satisfied  for  a  loaded  (i.e. 
dielectric)  body  of  revolution  (BOR) .  The  subroutine 
TESTSPHERE  generates  a  spherical  BOR  to  compare  numerical 
solutions  of  the  main  program  and  function  subprograms  with 
known  analytical  solutions.  TESTSPHERE  is  used  solely  for 
validating  the  subroutines  that  calculate  the  circular 
aperture  radiation  field.  It  is  not  used  in  the  calculation  of 
the  radome  scattered  fields.  In  all  configurations  the  antenna 
is  centered  on  the  origin  of  the  global  coordinate  system 
which  is  the  cylindrical  coordinate  system  of  the  BOR. 

B.  DEFINITION  OF  OGIVE 

The  ogive  is  a  figure  of  revolution  obtained  by  rotating 
an  arc  of  a  circle  about  an  axis  in  the  plane  of  the  arc  as 
depicted  in  Figure  2.1.  The  curve  of  the  ogive  is  rotated 
around  the  z  axis  to  produce  a  surface  of  revolution. 
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The  profile  of  the  surface  is  generated  in  the  r-z  plane 


by  the  equation 


r(z)  •jR2-zi*b-R 


(2.1) 


OGIVE  body  of  revolution  (bor),  r(z)  is  rotated 
around  z  axis  The  radius  (R)  of  tne  parent 
circle  originates  in  the  z  ■  0  plane. 


Figure  2.1  Ogive  as  Defined  by  Equation  (2.1). 


where  R  is  the  radius  of  the  parent  circle  and  b  is  the  radius 


of  the  base. 


C.  MODIFIED  OGIVE  FOR  LDBORMM.F 

In  the  program  LDBORMM.F  the  antenna  is  located  at  the 
origin  of  the  coordinate  system.  In  a  real  system,  however, 
the  antenna  may  be  displaced  on  the  z  axis,  from  the  origin  of 
the  defining  Equation  (2.1).  The  subroutine  OGIVE  prompts  the 
user  to  enter  the  value  z',  which  is  the  displacement  on  the 
z  axis  of  the  antenna  location  and  the  position  on  the  z  axis 
where  the  radius  (R)  of  the  parent  circle  of  the  ogive 
commences  rotation.  Figure  2.2  depicts  the  modified  profile  of 
an  ogive  as  calculated  by  the  subroutine  OGIVE. 


Figure  2.2  BOR  Generated  by  Equation  (2.2). 


The  modified  definition  for  the  radius  as  a  function  of 
position  on  the  z  axis  is  given  by, 
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r{z) 


b 


(2.2) 


zzz' 

Jr2- (z-z')2+b-R  z>z' 

where  z'  is  the  displacement  on  the  z  axis.  Using  this 
modified  definition  the  source  was  located  at  the  origin  of 
the  cylindrical  coordinate  system  used  to  generate  ogive(i.e., 
the  source  is  located  in  plane  z=0) .  To  reiterate,  the  z '  of 
Equation  (2.2)  and  Figure  2.2  denotes  the  distance  from  the 
origin  where  the  radius  of  the  parent  circle  is  defined  and  is 
not  related  to  the  source  coordinates. 

D.  SUBROUTINE  OGIVE 

The  subroutine  OGIVE  returns  values  for  the  global 
variables  NP,  ZH  and  RH.  NP  is  the  number  of  surface 
generating  points  on  the  ogive  in  the  r-z  plane.  In  the 
subroutine  the  range  of  NP  is  limited  to  3  <  NP  <  400.  The 
program  prompts  the  user  for  the  following  inputs:  radius  of 
the  parent  circle  (R) ,  z' ,  and  the  base  radius  (b) . 

A  standard  rule  of  thumb  for  the  convergence  of  the  series 
approximation  for  the  surface  current  is  that  the  spacing  of 
the  surface  generating  points  be  approximately  one-tenth 
wavelength  or  less.  The  following  algorithm  computes  the 
number  of  points  required  for  a  segment  size  of  approximately 
one-tenth  of  a  wavelength  on  the  BOR 
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Z  =  \J2  *R*b-bz 
0  =  ARCSIN  1 

(2.3) 

AL  =  B*R  +  Z' 

NP  =  2  *  INTEGER  (5  *AL)  +1  . 

NP  is  the  global  variable  for  the  number  of  points  in  the  r-z 
plane  of  the  ogive  [reference  Figure  2.2].  For  a  circular 
aperture,  the  field  incident  on  the  radome  tends  to  vary  more 
rapidly  as  the  angle  off  boresight  increases,  thus  more 
generating  points  are  needed.  By  slightly  modifying  the  above 
algorithm,  a  surface  with  any  arbitrary  segment  size  can  be 
generated. 

Equation  (2.3)  gives  the  algorithm  employed  by  the 
subroutine  OGIVE  to  calculate  the  number  of  points  (NP)  on  the 
ogive  surface  in  the  r-z  plane.  Data  generated  uniformly  along 
the  radome  surface  will  be  nonuniform  in  z,  r,  and  0.  The 
subroutine  TESTSPHERE  generates  data  points  with  a  uniform 
spacing  in  angle  0,  whereas  OGIVE  points  are  uniform  along  S. 
Thus  a  comparison  of  the  two  outputs  will  not  yield  coincident 
data  points.  Appendix  A  contains  the  source  code  for 
subroutines  OGIVE  and  TESTSPHERE.  Appendix  E  provides 
mathematical  details  relevant  to  the  distribution  of  data 
points  in  angle  0  for  the  subroutines  OGIVE  and  TESTSPHERE. 
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If  the  number  of  points  is  greater  than  400,  array  limits 
are  exceeded  and  the  main  program  will  write  a  warning  to  the 
screen  and  stop  the  program.  For  bodies  of  revolution  where 
the  arc  length  is  greater  than  40  \  the  dimensions  of  the 
following  global  variables  must  be  changed  accordingly 
throughout  the  program:  RH,  ZH,  Z,  R,  B,  C,  ZLO  and  ZL.  ZHn 
and  RHn  are  the  product  of  the  local  coordinates  zn,  rn 
(computed  by  the  BOR  subroutine)  and  the  wave  number  k=27r/A.. 
All  spatial  dimensions  throughout  the  main  program, 
subroutines  and  functions  are  in  wavelengths  (A). 

For  surface  shapes  with  definitive  geometries  other  than 
an  ogive  or  sphere  the  user  must  provide  an  appropriate 
subroutine.  The  subroutine  must  return  the  global  variables 
NP,  ZH,  RH,  b,  R  and  Z'.  These  variables  are  required  to 
define  the  generating  points  of  the  BOR  profile  (see  Eqn.2.2) . 
The  call  for  the  subroutine  must  be  inserted  in  the  logical 
block  of  code  in  the  main  program  which  establishes  the  BOR 
geometry.  The  subroutine  is  appended  to  the  main  program. 
Chapter  V  addresses  modification  of  the  main  program  for 
various  bodies  of  revolution.  The  program  LDBORMM . F  saves  the 
inputs  from  the  BOR  subroutine  and  writes  them  to  the 
external  sequential  file  OUTLDBOR . 


14 


III.  METHOD  OF  MOMENTS 


A.  OVERVIEW 

In  order  to  calculate  the  scattered  field  pattern  of  a 
dielectric  radome  the  electric  current  (Js)  must  be  integrated 
over  the  surface  of  the  body  of  revolution.  The  magnetic 
vector  potential  (A)  and  scalar  electric  potential  ($)  are 
defined  in  terms  of  the  surface  current 

£f(i?,e,4>)  =  -jo>4  (J^  -  =1^)  (3.1) 

where  Es(R,0,0)  is  the  scattered  electric  field  in  spherical 
coordinates,  and  L(»)  is  an  operator  introduced  for  notational 
convenience.  The  potentials  are 

dS>  (3*2) 

s 


*&r7faJf^ds'  <3.3) 

S 

0  =  — J-  Ve  .  (3.4) 

R  is  the  magnitude  of  the  difference  r  and  £',the  position 
vectors  of  the  field  and  source  points  respectively.  The 
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operator  ?,•  is  the  surface  divergence  on  S  [Ref.  6:  pp.  15- 
16].  Figure  3.1  illustrates  the  quantities. 


P(r,e,$> 


To  calculate  the  electric  field  ,  due  to  a  source 
with  an  electric  current  density  Jf 
Equations  (3. 1 )  through  (3  4)  must  be  calculated. 

Figure  3.1.  Scattering  due  to  Surface  Currents. 

The  objective  is  to  solve  for  the  current  on  the  surface 
of  the  radome.  Following  the  method  of  moments  (MM),  the 
electric  current  (J„)  on  the  surface  (S)  is  represented  by 

iZs  (*nj  Jai  +  laj  •  (3.5) 

in  Equation  (3.5)  and  are  the  tangential  and 

azimuthal  components  of  the  current  on  the  surface  of  the 
radome.  The  expansion  functions  for  are  chosen  to  be 


L6 


pulse  function 

Figure  3.2  Basis  Functions  Pulse  (P)  and  Triangle  (T)  of 
Equations  (3.6)  and  (3.7) 
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j‘  =  t^  eJa*  7  =  1,2, - WP-2  (3.6) 

ri 

e1a*  7=1,2, - iVP-1  (3.7) 

rj 

where  the  azimuthal  mode  number  is  n=0,±l,±2. . . .±«.  An  exact 
solution  requires  an  infinite  number  of  terms  in  the  sum  over 
the  index  n.  In  practice,  when  the  source  (antenna)  is  on  the 
axis  of  symmetry,  the  sum  converges  rapidly  and  only  a  few 
azimuthal  modes  are  required. 

The  unit  vectors  t  and  $  are  in  the  tangential  and 
azimuthal  directions.  The  functions  Tj(t)  and  Pj(t)  are  the 
triangle  and  pulse  functions  as  shown  in  Figure  3.2.  The 
abscissa  t  is  the  arc  length  along  the  generating  curve  of  the 
BOR.  It  is  assumed  that  the  generating  curve  consists  of  NP-1 
straight  line  segments  where  NP  is  an  integer.  The  expansion 
functions  of  Equations  (3.6)  and  (3.7)  are  especially 
appropriate  if  the  BOR  is  an  infinitely  thin  perfectly 

A 

conducting  surface  with  edges.  This  is  true  because  the  t 
directed  electric  current  approaches  zero  at  an  edge,  whereas 
the  0  directed  electric  current  may  grow  large  [Ref.  7:  p.  6]. 
The  MM  solution  of  the  EFIE  reduces  the  integral  equation  to 
a  matrix  equation  for  the  unknown  coefficients  I1  and  I*  .  This 
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is  accomplished  by  the  MM  testing  procedure.  The  test 
functions  (W^,^/)  are  chosen  to  be  the  complex  conjugates  of 
Equations  (3.6)  and  (3.7).  Equation  (3.5)  is  substituted 
into  Equations  (3.2)  to  (3.4)  for  the  solution  of  Equation 
(3.1).  The  dot  product  of  Equation  (3.1)  is  then  taken  with 
the  test  functions  These  dot  products  are  then 
integrated  over  the  surface  S.  The  resulting  matrix  equation 
is 


■j 

■z“ 

z# 

-1 

vc 

J* 

z* t 

z 

where  Z  is  a  square  matrix  called  the  impedance  matrix,  and  V 
is  a  column  vector  called  the  excitation  vector  [Ref.  7:  pp. 
6-31]  . 

The  impedance  elements  are  identical  to  those  described  by 
Mautz  and  Harrington.  Detailed  equations  for  them  appear  in 
reference  [Ref.  7]  and  will  not  be  repeated  here.  However,  the 
excitation  elements  will  depend  on  the  characteristics  of  the 
antenna  and  the  location  of  the  radome.  The  derivation  of  the 
excitation  elements  for  a  BOR  in  the  near  field  of  a  source 
has  not  been  presented  elsewhere,  and  therefore  will  be 
derived  in  detail  in  section  (3.B.). 

In  LDBORMM.F,  the  subroutine  ZMAT  calculates  the  impedance 
matrix  and  the  subroutine  GENEX  calculates  the  excitation 
vector  (V) .  The  subroutines  DECOMP  and  SOLVE  solve  the  system 
of  equations  resulting  in  the  electric  current  coefficients 
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(I*.,*)  .  The  current  coefficients  are  substituted  in  Equation 
(3.2),  which  determines  the  surface  current  (Js)  [Ref.  7: 
pp. 41-64 ]  . 

With  the  surface  current  (J.)  calculated,  the  scattered 
electric  field  (Es)  is  determined  by  Equation  (3.1).  The  total 
electric  field,  E(R,0,0),  is  then  given  by 

£(J?,6,4>)  =  El  *  El  <3-9> 

where  Ea  is  the  electric  field  of  the  circular  aperture. 

B.  NEAR  FIELD  OF  A  CIRCULAR  APERTURE 

In  order  to  maintain  a  clear  picture  of  the  analysis  that 
follows,  it  is  advantageous  to  keep  in  mind  the  following: 

1.  The  electric  current  distribution  over  the  circular 
aperture,  which  represents  an  ante^?  _ ,  is  defined  in  Cartesian 
coordinates . 

2.  The  electric  field  due  to  the  aperture  and  the  scattered 
electric  field  due  to  the  radoroe  are  computed  in  spherical 
coordinates  (R,0,<£). 

3.  Since  the  ogive  is  a  body  of  revolution  it  is  generated  in 
a  cylindrical  coordinate  system  (r,0,z). 

4.  Primed  quantities  denote  source  coordinates,  unprimed 
quantities  denote  coordinates  at  the  point  of  observation 
(with  the  exception  of  z',  which  was  discussed  earlier). 

To  calculate  the  excitation  vector  (V)  of  Equation  (3.8), 
the  electric  field  of  the  antenna  must  be  known.  Since  the 


20 


antenna  is  typically  operating  in  the  near  field  of  the 
radome,  far  field  approximations  are  not  valid  and  a  more 
accurate  method  of  calculation  must  be  employed.  Since  many  of 
the  radar  antennas  used  in  airborne  applications  are  circular 
paraboloids  or  slotted  waveguide  plates,  the  radiation  source 
is  modeled  as  a  circular  aperture  with  a  uniform  current 
distribution  polarized  in  the  x  direction 

'Kx'.y'.z')  =  JxX  .  (3.10) 

For  simplicity  the  electric  current  distribution  over  the 
aperture  is  assumed  to  have  constant  amplitude  and  a  linear 
phase 

Jx  =  joeJ**/sine«COB*.  .  (3.11) 

J0  is  the  amplitude  of  the  electric  current,  while  0S  and 
are  scan  angles  in  the  spherical  coordinate  system.  The 
coordinates  of  a  point  in  the  plane  of  the  aperture  are  given 
by  (x' ,y' ,0) . 

To  compute  the  Cartesian  components  of  the  electric  field, 
Jx  is  used  in  Equations  (6-108a)  to  (6-108c)  of  Reference 
[Ref.  5:  p.  283].  A  straightforward  transformation  of 
rectangular  to  polar  coordinates  (see  Appendix  E)  yields  the 
following  expressions  for  the  Cartesian  components  of  the 
scattered  electric  field 
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where 


E  =  ~Jr]  J° 
x  4nJc 


JJeXJ  [G1  +  (r  cos4>-x' cos4>/)2^2]r/  dr7  dtf  (3.12) 


y  4  7t  ic 


^^2  j J exl  (r  cosfy-z'  costy) 

(r  sin^-r7  sin^C^r7  dr'dfy' 


(3.13) 


_  -J'n^o 


j J  exl(r  cosfy-r'  cosfy)  z  G2  r'  dr'  d<J>7  (3.14) 


xl  =  jk  (r7  cos<t>7sin0a-l?) 


(3.15) 


i?  =  J  (r  -r7)2  +  z2  +4  rr 


(3.16) 


Gi  =  k2R2-l-jkR 
R3 


(3.17) 


3  +3jkR-k2R2 
R5 


(3.18) 


Figure  3.3  defines  the  quantities  used  in  the  equations.  It  is 
assumed  that  the  antenna  is  scanned  such  that  <ps=  0. 

The  rectangular  components  of  the  electric  field  are 
converted  to  spherical  components  using  the  transformation 
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Eq  ■  Excos8cos^  +  £ycos0sin^-£tsin0 


(3.19) 


Ek  -  Exsindcos$  +  Eysin&9in$  +  EacosQ 


E+  -  -£xsin$ +  £yco8$  . 


The  function  subprograms 

CIRCTHETA  ,  CIRCRHO  and  CIRC  Phi  calculate  the  spherical 
electric  field  components  at  the  point  (  P  ) 

The  point  (P)  Is  defined  in  cylindrical  coordinates, 
r(z),  0  and  Z 

The  primed  coordinates  are  In  the  plane  Z  •  0 
The  circular  aperture  Is  centered  at  the  origin 

Figure  3.3  Geometry  for  Equations  (3.12)  to  (3.18) 

The  function  subprograms  CIRCTHETA,  CIRCRHO,  and 
calculate  the  spherical  electric  field  components  E#  , 
E#,  respectively. 


(3.20) 

(3.21) 


CIRCPHI 
Er  ,  and 


1.  FUNCTION  SUBPROGRAMS,  CIRCTHBTA,  CIRCRHO,  AND  CIRCPHI 

The  function  subprograms  CIRCTHETA,  CIRCRHO  and 
CIRCPHI  are  called  in  the  subroutine  GENEX.  In  the  main 
program,  the  vector  field  due  to  the  source  (Ea)  is  added  to 
the  scattered  field  vector  (Es)  of  the  body  of  revolution  to 
obtain  the  total  field  (ER 6  ,  per  Equation  (3.9).  Since  the 
far  field  approximation  is  valid  in  this  case,  the  closed  form 
solution  for  a  circular  aperture  derived  in  Section  III.B.l.b 
can  be  used.  Subroutine  GENEX  requires  the  components  of  Ea  in 
the  near  field.  The  function  subprograms  calculate  the 
spherical  electric  field  components  due  to  a  uniformly 
illuminated  circular  aperture  at  an  arbitrary  point  in  space 
specified  by  the  spherical  coordinates  R,  0,  0.  The  source  can 
be  scanned  off  boresight  by  entering  non-zero  scan  angles  when 
the  program  is  executed. 

In  the  subroutine  GENEX,  these  functions  are  called 
for  each  point  of  integration  in  the  tangential  and  azimuthal 
components  of  the  excitation  vector  (V)  which  occur  in 
Equation  (3.8),  and  are  described  in  the  following  section. 

Appendix  D  contains  a  table  of  the  variables  in  the 
argument  list  of  the  function  subprograms  and  the  source  code 
for  CIRCSUB . F .  The  program  CIRCSUB.F  was  used  to  validate 
the  algorithms  of  CIRCTHETA,  CIRCRHO  and  CIRCPHI.  The 
programs  TESTCIRC.M  and  SCANPLT.M  (Appendix  D)  were  used  to 
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generate  the  analytic  solution  and  plot  the  comparative 
results. 

a.  Numerical  Method 

To  determine  the  electric  field  of  the  antenna  and 
the  excitation  vector  requires  evaluating  several  integrals. 
In  general  these  integrals  can  not  be  reduced  to  a  closed  form 
expression  for  the  geometries  under  consideration.  The 
integrals  must  therefore  be  integrated  numerically.  The  chosen 
method  for  evaluating  all  integrals  in  this  program  is 
Gaussian  quadrature.  In  Gaussian  quadrature  algorithms,  the 
integral  is  over  a  normalized  interval,  such  that  -1<  x  <  1. 
N-point  Gaussian  quadrature  approximates  an  integral  by 

f,  f(X)  <*X  *  EL  <■>»  f(X„)  >3.22) 

where 


N  =  2/n  +  2  (3.23) 

and  CiJ^  are  the  coefficients  and  Xn  are  the  zeros  of  the 
Legendre  polynomial  of  order  m  (m=0,+l,±2 . .  .<»)  [Ref.  8]. 

The  number  of  terms  in  the  sum  is  always  an  even 
number.  The  program  GAUS.F  (see  Appendix  C.)  will  generate  an 
external  sequential  file  for  any  range  of  N  where  0  <  N  <  500. 
The  first  record  in  the  file  is  an  even  number  (N)  .  The 
remainder  of  the  file  contains  N  rows  of  xn  and  Hi  •  The  gaus 
files  are  stored  in  a  subdirectory  named  "gaus"  of  the 
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directory  containing  the  main  program  LDBORMM.F.  The  path  to 
the  data  files  is  "(MAIN  DIRECTORY)/gaus/gaus### ,  where  ###  is 
N.  To  create  the  data  file  (gaus###) ,  execute  GAUS.F  in  the 
subdirectory  "gaus."  The  screen  prompts  the  user  for  the  data 
filename  and  an  even  number  of  points. 

For  integrals  that  are  over  a  range  other  than  -1 
to  +1,  a  change  of  variables  is  necessary  before  Gaussian 
quadrature  is  applied.  Integrations  over  the  antenna  aperture 
(from  zero  radius  to  the  antenna  radius  rD)  require  this 
change  of  variables  is  given  by  [Ref  8:  pp.  523-527] 


A  similar  technique  is  used  to  evaluate  the  <t> 
integrals.  The  interval  of  integration  for  0  is  symmetric  and 
does  not  require  a  change  of  variables  since,  -ir<  <p  <  n 

^n-KX„*%  .  (3.25) 

Equations  (3.12)  through  (3.14)  are  integrated  by  this  method. 
In  the  function  subprograms,  CIRCTHETA,  CIRCRHO  and  CIRCPHI, 
the  variables  of  integration  over  the  aperture  are  PHIPRIME 
(0')  and  RP  (r')  and  the  limits  are  from  0  to  2 ir  and  0  to  r0, 
respectively. 
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In  order  to  verify  the  algorithms  and  coding  for 
CIRCTHETA,  CIRCRHO  and  CIRCPHI  the  field  magnitude  Ee  was 
computed  and  compared  with  the  analytical  solution  given  by 


j kra2 E0e~jk*  Jx(icr0  (sin0-sin05) 
“  R  krQ(  sine-sin8,) 


COS ( 0 )  (3.26) 


where  E0=r?Jo  is  the  electric  field  amplitude,  is  the  first 
order  Bessel  function,  0,  is  the  scan  angle  and  (R)  is  the 
distance  from  the  center  of  the  aperture  to  the  far  field 
point  (P)  of  Figure  3.3.  It  is  assumed  that  R  is  much  greater 
than  r'  in  the  above  formula.  Therefore  any  vector  from  a 
source  point  to  the  observation  point  is  approximately 
parallel  to  the  vector  R  [Ref.  9:  pp.  478-487]. 

The  far  field  distance  for  which  Equation  (3.26)  is 
valid  depends  on  radius  (r0)  of  the  aperture.  Antenna 
engineers  commonly  consider  the  far  field  to  exist  at 
distances  greater  than 


However  in  optics  the  far  zone  applies  to  values  of  R  where 
the  Fresnel  number  is  less  than  one.  Using  the  latter 
criterion  the  minimum  distance  for  which  Equation  (3.26)  is 
valid  is  given  by  [Ref.  10:  pp.  95-96] 
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(3.28) 


1  i 


kron> 


I”. 


2nR 

2 


The  program  TESTCIRC.M  in  Appendix  D  is  used  to 
calculate  the  analytical  solution  for  E9  in  accordance  with 
Equation  (3.26).  The  results  are  compared  with  the  calculated 
field  computed  by  CIRCSUB.F.  In  Figures  (3.4),  (3.6),  and 
(3.8)  the  analytical  solution  is  plotted  with  a  solid  line 
while  the  values  calculated  by  the  test  program  CIRCSUB.F 
are  plotted  as  points. 

Figures  (3.4),  (3.6),  and  (3.8)  show  the 
convergence  of  the  numerical  solution  from  CIRCSUB.F  as  a 
function  of  the  number  of  integration  points  with  the 
analytical  solution  given  by  Equation  (3.26)  for  one,  two  and 
five  wavelength  radius  circular  apertures.  Figures  (3.5), 
(3.7),  and  (3.9)  show  plots  of  the  field  magnitudes  (E9,  ER, 
Et)  versus  the  angle  in  degrees.  Figures  3.10  through  3.12 
show  plots  from  CIRCTHETA  and  Equation  (3.26)  for  various  scan 
angles  and  aperture  radii.  (Note  that  the  vertical  scale  in 
Figures  3.4  to  3.12  is  in  volts/meter,  not  dB.) 

Table  3.1  summarizes  parameters  for  Figures  (3.4) 
through  (3.12).  The  plotted  results  verify  the  excellent 
agreement  between  the  algorithm  of  the  function  subprograms 
and  the  analytical  solution.  It  is  important  to  note  the 
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number  of  points  of  integration  in  r '  (CNRHO)  and  <p'  (CNPHI) 
of  Equations  (3.13)  through  Equation  (3.15)  required  for 
convergence  to  the  analytical  solution  of  Equation  (3.26). 
The  number  of  points  of  integration  in  the  table  are  the 
minimum  number  required  for  convergence.  As  the  radius  of  the 
aperture  increases,  the  number  of  iterations  required 
increases  as  the  product  of  CNRHO  and  CNPHI.  For  example, 
Test  3  requires  1200  iterations  for  convergence  to  the 
analytic  solution.  In  the  program  LDBORMM.F  the  functions 


TABLE  3.1.  DATA  SUMMARY  FOR  TEST  CASES. 


Test  Radius 
.  (r0) 


subprograms  are  called  five  times  in  the  subroutine  GENEX 
which  is  in  turn  called  once  for  each  azimuthal  mode.  Thus  the 
number  of  calls  for  the  function  subprograms  is  the  product  of 
the  number  of  integration  points  in  t  (NT) ,  number  of 
integration  points  in  <p  (NPHI)  ,  the  number  of  azimuthal  modes 


(MODES)  and  the  number  of  generating  points  on  the  surface  of 
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the  BOR.  For  a  large  BOR,  the  computer  run  time  can  become 
significant.  It  is  therefore  prudent  to  test  the  function 
subprograms  for  the  minimum  number  of  points  of  integration 
required  for  convergence  before  executing  the  main  program. 
2.  THE  SUBROUTINE  GENEX 

The  subroutine  GENEX  calculates  the  excitation  vector 
(V)  for  Equation  (3.8)  .  All  variables  in  the  argument  list  are 
inputs  with  the  exception  of  the  excitation  vector.  Appendix 
B  contains  a  list  of  variables  in  the  argument  of  the 
subroutine  GENEX. 

The  jth  element  of  V1  and  V*  are  given  by 

V^p  =  -  f  f  Wfin  -  El  ds,  j-1,2,  .  .  . NP-2  (3.29) 

■n  J  g  — 

±  ffswtn- El  ds,  j-1,2 - NP-1  (3.30) 

where  the  indices  j  and  n  are  the  same  as  Equations  (3.6)  and 

(3.7) ,  and  WL*  is  the  testing  function.  The  test  functions  are 
chosen  by  Galerkin's  method  and  therefore  are  the  complex 
conjugates  of  the  expansion  functions  Q}4  (Equations  (3.6)  and 

(3.7) ].  Epa  is  the  pth  component  electric  field  of  the 
antenna  and  S  denotes  the  surface  of  the  radome.  The 
representation  p  is  a  generalized  for  the  spherical  components 
R,  0  and  <p  .  These  components  are  provided  by  the  subroutines 
CIRCTHETA,  CIRCRHO,  and  CIRCPHI . 
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Figure  3.4  Test  1. 
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Figure  3.6  Test  2. 
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Figure  3.7  Test  2,  E(R,0,«)  Calculated  by  CIRCSUB.F. 
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field  magnitude  vs. theta 


Figure  3.9  Test  3,  E(R,0,0)  calculated  by  CIRCSUB.F. 
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TEST  OF  SUBROUTINE  CIRCTHETA 


Figure  3.10  Test  4,  0  Scanned  to  30  Degrees. 
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TEST  OF  SUBROUTINE  CIRCTHETA 
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Figure  3.11  Test  5,  0  Scanned  to  45  Degrees. 
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TEST  OF  SUBROUTINE  CIRCTHETA 


Figure  3.12  Test  6,  6  Scanned  to  60  Degrees. 
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The  quantities  ER,  Es  and  E#  completely  define  £‘  and  are  used 
by  GENEX  to  evaluate  Equations  (3.28)  and  (3.29)  numerically. 
To  further  reduce  equations  (3.28)  and  (3.29)  to  a  form 
suitable  for  programming  requires  evaluation  of  the  dot 
products.  Referring  to  Figure  3.13,  the  tangent  direction 
along  the  surface  can  be  written  as 

£  =  f  sinVj  +  £  cosvj  (3.31) 

By  the  orthogonality  properties 

$  •  £  *  0 

|  • 0  =  0  (3.32) 

t  •  $  *  0  . 

Therefore,  from  (3.28)  and  (3.29)  vt¥=vw=v^=0.  The  remaining 
elements  are 


Vnt=  f  f  T*  CFx  Er  dt  +  f  TJ  CFU1  Er  at  Id*  (3.33) 

J  4  [•'A, 

Vn?  =  f  f  T*  SF  i  Eq  dc  *  f  T]  SF ^  (3.34) 

v*f=[  e -jn*  f  Pi  dt  d<t>  (3.35) 

A,  =  t;.,  -  t;  (3.36) 
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Ai„i  =  ti*  2  ~  1 


(3.37) 


SFt  =  sin (Vj-  00  (3.39) 

CFX  =  COS  (Vi-  0')  •  (3.40) 

The  quantities  are  defined  in  Figures  3.2  and  3.13.  The 


Figure  3.13  Geometry  for  Equations  (3.30)  through 
(3.39) 


elements  of  the  excitation  vector  (V)  are  stored  in  the  global 
variable  R,  where 
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VitS+VitR  is  in  R(i)  i=l  to  NP-2 

V **  is  in  R ( i+MT)  i=l  to  NP-l 

and  MT=NP-2  [Ref.  7:  p.  57]. 

C.  MODIFICATION  OF  THE  IMPEDANCE  MATRIX  FOR  DIELECTRICS 
1.  THIN-SHELL  APPROXIMATIONS  FOR  DIELECTRICS 

For  a  perfect  conductor  the  surface  resistance  (Rs)  is  zero. 
To  satisfy  the  boundary  condition  requires 

^  =  -E*  (3.41) 

for  the  tangential  components  of  the  incident  and  scattered 
fields.  In  this  case  the  incident  field  is  that  due  to  the 
antenna.  When  a  dielectric  body  is  present,  the  incident 
electric  field  produces  a  polarization  current  rather  than  a 
conduction  current.  In  this  case  the  boundary  conditions  state 
that  the  electric  and  magnetic  fields  must  be  continuous.  The 
total  electric  and  magnetic  fields  are 


E  =  El  +  El 
U  =  01  +  HI  . 

Maxwell's  equations  must  be  satisfied  where 

V  x  £  =  -ju\i0  H 


(3.42) 


(3.43) 


and  by  superposition 
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v  X  El  -  -JU|i0  HI 


(3.44) 


so  that 


V  X  H  =  JO)C0  £  +  ju  (c  -e0)  £  (3.45) 

where  the  second  term  on  the  right  side  of  Equation  (3.44)  is 
nonzero  for  t*tQ  . 

The  polarization  current  may  be  defined  as 

Jg  s  jo)  (e-e0)  £  .  (3.46) 

Maxwells'  second  equation  becomes 

V  X  tfg  =  ju>t0E^  *  Jp  .  (3.47) 


Thus  Es  and  Hs  can  be  determined  from  the  polarization  current 
(Jp)  radiating  in  free  space.  For  dielectric  bodies  with  er>>l 
or  very  thin  dielectric  shells  with  thickness  t  <<  A,  the 
tangential  component  of  the  polarization  current  will  be  much 
greater  than  the  normal  component.  In  these  cases  the  normal 
component  can  be  neglected  and  a  modified  form  of  the  EFIE 
results 


L(JP)  +■ 


El 


(3.48) 


JO)  (e  -e0)  t  - 

where  the  coefficient  of  Jp  is  called  the  surface  loading  or 
surface  impedance  (Rs) .  (Frequently  it  is  referred  to  as  a 
surface  impedance  because  it  can  be  complex.  Here  it  will  be 
called  a  resistance  to  avoid  confusion  with  the  surface 
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impedance  approximation  in  which  both  electric  and  magnetic 
currents  are  allowed.)  For  a  lossless  material  Rs  is 
imaginary;  for  a  resistive  material  it  is  real. 

The  new  matrix  equation  corresponding  to  Equation 
(3.8)  is 


[J]  =  IZ"  +  ZL  )-'  [V\  (3.49) 

where 


Zm  =  ffsL ta'UJj  &  (3.50) 

=  ffjLn  ‘  RsJ^ds  .  (3.51) 


2.  EVALUATING  THE  IMPEDANCE  MATRIX 

The  impedance  matrix  (Z)  of  Equation  (3.8)  is 
calculated  by  the  subroutines  ZMAT,  ZLOAD,  and  ZTOT.  If  for 
some  reason  the  impedance  of  the  radome  is  zero  (a  perfect 
conductor) ,  only  the  original  impedance  matrix  from  ZMAT  is 
needed.  For  other  materials  the  load  impedance  is  calculated 
by  ZLOAD  and  added  to  the  ZMAT  result  the  subroutine  ZTOT. 

The  thin  shell  approximation  assumes  the  thickness (t) 
is  much  less  than  the  wavelength  (A)  and  the  skin  depth  (6)  of 
the  shell  material.  This  is  a  valid  assumption  for  most 
practical  radome  designs  and  materials.  The  polarization 
current  is  primarily  tangential  to  the  surface  of  the 
dielectric  and  the  normal  components  are  negligible.  This 
approximation  is  applicable  to  lossy  as  well  as  loss  free 
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dielectrics.  In  the  lossy  case  the  load  matrix  (ZL)  has  a  real 
and  imaginary  component,  while  the  loss  free  case  is  purely 
reactive. 

At  this  point  it  is  necessary  to  relate  the  surface 
impedance  to  the  electrical  properties  of  the  material  such  as 
the  dielectric  constant  and  loss  tangent.  The  electric  field 
due  to  a  thin  dielectric  shell  is  given  by  equation  (3.48) 
with  Rs  defined  by 


r  =  _ — ± - 

jo)  A  e  t 


(3.52) 


where  Ae=e-e0  [Ref.  11:  pp.  531,532].  As  stated  earlier,  the 
dielectric  constant  of  the  shell  may  be  complex.  Letting 
e=e'-je"  gives 


*5  = 


where  e  '  =en  e. 


jo)  t 

The  electric  loss  tangent  is  given  by 


(3.53) 


tan  (6)  =  —  =  —5—  (3.54) 
t' 

where  <5  is  the  loss  angle  and  a= ue"  is  the  conductivity 
representing  all  losses  in  the  medium.  Therefore  the  loss 
tangent  is  a  measure  of  the  power  loss  in  the  medium.  A  medium 
is  a  good  conductor  if  a» coe ,  and  a  good  insulator  if  a« <oe 
[Ref.  2:  pp.  342,343].  For  <oe0  ■  1/60 A,  and  thickness  t=nA 
substitution  of  (3.58)  into  (3.53)  yields 


45 


(3.55) 


s  n  [  j  (er-l)  +ertanfi] 

Typical  radome  materials  have  2  <  er  <  10  and  0.0005  <  tan  <5 
<  0.005  .  Inhomogeneous  radomes  (i.e.  sandwiches)  can  be 
treated  by  using  an  equivalent  dielectric  constant. 
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IV.  PROGRAM  DESCRIPTION 


A.  DATA  FLOW 

The  program  LDBORMM.F  has  eight  subroutines  and  four 
function  subprograms.  Subroutines  and  functions  not  covered 
here  are  documented  in  reference  [Ref.  7].  Table  4.1  gives  a 
brief  description  of  all  subroutines  and  function  subprograms 
in  LDBORMM.F. 


TABLE  4.1  SUMMARY  OF  SUBPROGRAMS 


NAME 

TYPE 

DESCRIPTION 

•••  : 

OGIVE 

subroutine 

defines  the  radoxne  dimensions 

GENEX 

subroutine 

calculates  the  excitation  vector 

ZMAT 

subroutine 

calculates  the  moment  matrix 

ZLOAD 

subroutine 

calculates  the  submatrix  for  the 
dielectric  BOR 

ZTOT 

subroutine 

calculates  the  moment  matrix  for 
dielectric  BOR 

DECOMP 

subroutine 

inverts  the  moment  matrix 

SOLVE 

subroutine 

calculates  the  current 
coefficients 

PLANE 

subroutine 

calculate  E  far  field 

BLOG 

function 

fourth  order  series  expansion  of 
log(x) ,  for  x<  0.1. 

CIRCTHETA 

CIRCRHO 

CIRCPHI 

functions 

calculate  the  p,  0  and  <p 
components  of  the  electric  field 
at  a  point  (P) . 

The  block  diagram  of  Figure  4.1  shows  the  structure  of  the 


calls  for  subroutines  and  functions  in  the  main  program. 


Figure  4.1  Structure  of  main  program  LDBORMM . F . 
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B.  PROGRAM  PARAMETERS 


1.  NUMBER  OF  POINTS  ON  BOR 

The  number  of  points  (NP)  must  be  at  least  3  but  less 
than  or  equal  to  400.  This  restriction  is  due  to  the  array 
dimensions  set  in  the  code  and  is  not  a  limitation  in  the 
solution.  The  subroutine  OGIVE  generates  segments  on  the 
surface  of  the  BOR  at  approximately  one  tenth  of  one 
wavelength.  Thus,  a  BOR  with  an  arc  length  greater  than  40  X, 
matrix  dimensions  are  exceeded  for  the  global  variables  RH, 
ZH,  Z,  R,  B,  C,  ZLO,  and  ZL.  The  subroutine  TESTSPHERE  allows 
the  user  to  select  the  number  of  points  on  the  surface  of  a 
sphere.  The  points  generated  by  TESTSPHERE  are  spaced 
uniformly  in  the  angle  0. 

2 .  DIRECTORY  STRUCTURE 

a .  Subdirectory  GAUS 

A  subdirectory  named  "gaus"  must  contai  l  the 
external  sequential  files  with  the  weights  and  abscissas  for 
the  Gaussian  quadrature  algorithm.  The  data  files  are 
generated  by  the  program  GAUS.F,  which  is  executed  in  the 
subdirectory  "gaus."  Appendix  C  contains  the  program  GAUS.F. 
The  status  of  the  files  in  the  main  program  is  "OLD",  thereby 
invoking  an  error  condition  if  the  file  does  not  exist.  Figure 
4.2  illustrates  the  directory  structure. 
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External  Input/Output 


The  external  sequential  data  files  required  for  the 
main  program  are  shown  in  Figure  4.2.  The  programs  that 
generate  the  files  are  indicated  by  dashed  lines  and  arrows. 
The  dashed  line  originates  at  the  generating  program  and  the 
arrow  terminates  on  the  file  generated. 

c.  The  Program  TEST.M 

The  programs  TEST.M  and  TESTCIRC.M  must  be  executed 
before  the  plot  routines  GENPLT.M  and  SCANPLT.M.  The  program 
TEST.M  (or  TESTCIRC.M)  generates  the  analytical  solution  given 
by  Equation  (3.26).  The  parameters  for  BOR  radius,  antenna 
radius  and  scan  angle  are  entered  by  the  user  .  These 
parameters  must  be  the  same  as  the  inputs  to  the  main  program. 
The  routine  GENPLT.M  plots  the  analytical  solution  generated 
by  TEST.M  and  the  numerical  solution  calculated  by  the  main 
program  for  comparison.  A  comparison  of  the  two  results 
illustrates  the  effect  of  the  radome  on  the  antenna  pattern. 

d.  The  Subdirectory  CIRCTEST 

The  subdirectory  CIRCTEST  contains  the  programs  and 
data  files  required  to  test  for  the  number  of  terms  necessary 
for  the  series  approximation  of  the  EFIE  to  converge  to  the 
analytic  solution.  The  main  program  requires  extensive 
execution  time  for  large  BOR  or  large  antenna.  In  order  to 
minimize  run  time,  the  minimum  number  of  integration  points 
over  the  antenna  is  determined  by  the  execution  of  the 


programs  in  the  subdirectory  CIRCTEST.  The  plotting  routine 
SCANPLT . M  plots  the  results  for  comparison. 

Table  4.2  provides  a  summary  of  the  files  in  this  thesis. 


TABLE  4.2  FILE 

SUMMARY 

ess  ia 

ldbormm. f 

main  dir./ 

main  program,  solves  for 
the  radiation  pattern  of 
radome  and  antenna 

outldbor 

main  dir./ 

external  data  file  for 
output  Of  LDBORMM. F 

bcpole.m 
bxpole .m 
bang.m 

main  dir./ 

data  generated  by 

LDBORMM. F  for  plotting 
with  GENPLT. M 

test .m 

main  dir./ 

program  to  generate 
solution  for  Eqn.(3.26) 

ctest .mat 

main  dir./ 

data  generated  by  TEST.M 

genplt .m 

main  dir./ 

plotting  routine 

gaus .  f 

main  dir./ 
gaus/ 

generates  gaus###  files 

gaus### 

main  dir./ 
gaus/ 

~ 

data  files  for  Gaussian 
quadrature  integration 

circsub. f 

main  dir./ 
circtest/ 

program  to  test  function 
subprograms  for  minimum 
number  of  integration 
points 

testcirc.m 

main  dir./ 
circtest 

generates  -  ration  for 

E c~ ■  '  ‘  . 

circ.mat 

main  dir./  .rated  by 

circtest  ,.:.M 

scanplt .m 

main  dir./ 
circtest 

plotting  routine 

C.  MAIM  PROGRAM  MODIFICATIONS 

The  main  program  is  designed  in  a  modular  fashion  in  order 
to  facilitate  modifications  for  other  radome  shapes,  antenna 
types  or  dielectric  profiles. 
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1.  Modifications  for  Radome  Shapes 

The  main  program  contains  subprograms  to  generate 
ogive-shaped  or  spherical  bodies  of  revolution.  When  executed 
the  main  program  prompts  the  user  for  the  type  of  BOR  desired. 
Additional  types  of  BORs  may  easily  be  added  to  the  existing 
menu.  The  following  modifications  are  required: 

(a)  Insert  code  in  main  program  to  write  selection  number  and 
type  of  shape  to  screen  menu. 

(b)  Insert  code  in  logical  block  that  calls  subroutine 
generating  selected  shape. 

(c)  Append  subroutine  program  to  main  program.  Subroutine  must 
return  the  variables  NP,  RH,  ZH,  b,  R  and  Z' . 

These  variables  are  required  for  subsequent  calculations  and 
data  entry  to  external  files. 

2.  Modifications  for  Antenna  Types 

To  modify  the  main  program  for  an  antenna  type  that  is 
not  a  circular  aperture: 

(a)  Substitute  the  alternate  function  name  in  the  main  program 
where  ETF  and  EPF  are  assigned.  These  are  the  far  field  0  and 
0  components  of  the  antenna,  so  a  closed  form  approximation 
can  be  used  if  one  exists. 

(b)  Substitute  the  alternate  functions  in  the  subroutine 
GENEX,  to  calculate  the  variables  SI,  S2,  S3,  S4  and  S5. 

(c)  Append  the  functions  to  the  main  program. 
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To  be  consistent  with  the  model  developed  earlier,  the 
function  subprograms  must  calculate  £(R,0,0)  in  spherical 
components  at  points  in  the  near  field. 

3.  Modification  of  Dielectric  Profile 

The  vector  ZLO  has  a  number  of  elements  equal  to  the 
number  of  segments  (NP-1)  on  the  surface  of  the  loaded  body  of 
revolution.  The  elements  in  the  array  are  the  complex  values 
of  the  surface  resistance  (RJ  on  each  segment  of  the  BOR.  As 
written,  the  main  program  prompts  the  user  to  enter  the 
complex  surface  impedance  when  the  main  program  is  executed. 
As  presently  coded,  the  entered  value  is  assigned  to  each 
segment.  In  order  to  examine  the  effects  of  a  nonuniform 
surface  impedance: 

(a)  Modify  the  program  to  read  an  external  sequential  file  of 
length  (NP-1) ,  with  the  values  of  the  desired  impedance 
profile . 

(b)  Store  the  values  in  the  array  ZLO  in  the  main  program 
before  calls  to  ZLOAD  and  ZTOT.  The  order  of  storage  in  ZLO 
must  correspond  to  the  order  of  the  coordinates  stored  in  RH 
and  ZH. 
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V.  SUMMARY  OF  RE8ULTS 


A.  VALIDATION  OF  PROGRAM 

In  order  to  validate  and  test  the  main  program,  the 
surface  resistance  was  assigned  a  very  large  value.  For  very 
large  values  of  surface  resistance,  the  scattered  electric 
field  (Es)  approaches  zero.  The  radiation  pattern  for  a  BOR  in 
the  far  field  with  very  large  surface  impedance  is  therefore 
the  far  field  pattern  of  the  antenna.  For  a  circular  aperture, 
the  far  field  radiation  pattern  is  given  by  Equation  (3.26). 

Figures  5.1  through  5.3  compare  the  values  for  the 
radiation  pattern  calculated  by  the  main  program  and  the 
analytic  solution  generated  by  TEST.M.  The  calculated  values 
are  plotted  as  (.)  and  the  closed  form  solution  to  (3.26)  is 
a  solid  line.  The  ordinate  axis  units  are  volts/meter.  The 
abscissa  is  the  angle  from  the  axis  of  symmetry  0  in  degrees. 
The  plots  of  Figures  5.1  to  5.3  show  the  magnitude  of  the 
electric  field.  The  actual  field  values  are  complex  and  it  is 
the  complex  values  that  are  used  to  calculate  the  excitation 
elements.  The  figures  show  the  parameters  NP,  NT,  NPHI,  CNRHO , 
CNPHI ,  SCAN,  and  ARAD  as  defined  in  Appendix  B.  The  angle  <p 
of  the  receiver  is  zero  for  all  test  cases.  The  number  of 
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points  of  integration  for  the  antenna  were  determined  by  the 
program  JIRCSUB.F  to  insure  that  aperture  integrations 
converged.  Figures  3.4,  3.6,  and  3.8  show  the  results.  The 
number  of  points  of  integration  needed  over  the  BOR  surface 
(NT  and  NPHI )  becomes  very  large  as  the  radial  extent  (hence 
circumference)  of  the  radome  increases.  In  addition  to  an 
increase  in  the  number  of  integration  points  required,  the 
number  of  azimuthal  modes  also  increases  for  a  large  BOR.  In 
all  test  cases  the  mode  number  is  one  (n=±l) .  Test  results 
indicate  more  points  of  integration  and  higher  modes  are 
required  for  the  calculated  values  to  converge  to  the  analytic 
solution  as  the  antenna  radius  or  sphere  radius  is  increased. 
The  execution  times  were  in  the  range  of  2  to  4  hours  on  a  Sun 
Sparcstation  2. 

Figures  5.4  through  5.6  show  plots  of  the  magnitude  of  the 
spherical  electric  field  components  for  the  three  test  cases. 
As  expected  the  scattered  field  components  ETSCAT  and  EPSCAT 
are  40  db  below  the  source  field  (ETF)  since  in  all  cases  Rs 
is  10000.  EPF  is  zero  since  the  antenna  radiation  is  0 
polarized  in  this  plane.  Figures  5.4  through  5.6  are  plots  of 
field  magnitude  versus  angle  0. 

The  plots  show  the  convergence  of  the  calculated  values  as 
measured  against  the  known  analytic  solution  of  a  circular 
aperture  in  the  far  field.  These  results  verify  the  program 
for  the  case  where  the  radome  is  nearly  transparent  for  large 
real  surface  resistance  in  accordance  with  Equation  (1.1). 
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Figure  5.4  Electric  Field  Components  Calculated  by  LDBORMM.F 
for  Test  1. 
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When  the  radome  is  transparent  the  calculated  radiation 
pattern  becomes  the  closed  form  solution  of  the  far  field. 

B.  EFFECTS  OF  RADOME  ON  TRANSMISSION  OF  ELECTRIC  FIELD 

Test  cases  were  conducted  for  various  radome  shapes  and 
surface  impedances.  In  Figures  5.7  through  5.9  the  calculated 
radiation  patterns  with  the  radome  are  compared  to  the 
patterns  of  the  isolated  antenna.  A  comparison  of  the  two 
curves  illustrates  the  effect  of  the  radome  on  the  antenna's 
performance.  The  following  parameters  were  constant  for  the 
three  test  cases: 

*  resistance  =  0+jl700  n;  (reflection  coefficient  r~0.2.). 

*  number  of  azimuthal  modes  =  2  . 

*  number  of  integration  points  ;NT=2  , NPHI=40 ,  CNRHO=20 ,  CNPHI=4 0  . 

*  antenna  radius  r0  =  2  A. . 

*  radius  of  the  BOR  base  was  3A  (open  base-no  body  structure)  . 
The  radome  shapes  tested  were: 

*  TEST  1.  a  sphere  of  radius,  R=3A.. 

*  TEST  2.  a  cone  with  che  half  angle  =  30°. 

*  TEST  3.  an  ogive  with  a  parent  circle  radius  R  =  7.5^. 

The  forward  direction  is  the  region  0  <  90°  while  the  rear 
hemisphere  is  0  >  90°.  The  plotted  data  in  Figures  5.8  and  5.9 
have  significant  backscatter  (0>9O°)  due  to  the  shape  of  the 
radome.  The  electric  field  magnitude  of  the  sphere  (Figure 
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5.7)  for  0>9O°  was  below  a  -20  db  threshold.  This  threshold 
was  arbitr  rily  established  for  purposes  of  plotting  the  data. 
The  backscatter  from  the  sphere  is  negligible  since  the 
radiation  from  the  antenna  at  each  point  on  the  surface  is 
near  normal  incidence.  The  cone  and  the  ogive  shapes  reflect 
energy  towards  other  parts  of  the  surface.  These  multiple 
reflections  raise  the  sidelobes,  especially  at  wide  angles 
from  the  antenna  main  beam. 

The  coupling  of  the  reflected  energy  to  the  antenna  is 
assumed  to  be  negligible  in  this  model.  This  assumption  is  not 
valid  for  geometries  or  dielectric  materials  which  cause 
significant  backscatter,  because  the  reflections  from  the 
radome  perturb  the  antenna  excitation.  Not  only  does  this 
effect  increase  the  sidelobes,  but  also  causes  an  increase  in 
the  input  VSWR.  For  a  well  designed  radome,  this  mutual 
interaction  is  a  secondary  effect. 

C.  RECOMMENDATIONS  FOR  FURTHER  STUDY 

The  following  suggestions  are  made  for  further  development 
and  improvement  of  the  program  LDBORMM.F  : 

1.  Obtain  test  data  for  a  real  system.  Run  the  program  for  the 
data  and  compare  the  results  of  the  computer  model  to  the 
actual  measured  test  data. 

2.  Develop  test  data  from  an  experimental  model  constructed  in 
the  ECE  Departmental  Transient  Electromagnetic  Scattering  Lab. 
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Compare  the  measured  fields  with  the  predictions  of  the 
computer  model. 

3 .  Perform  a  parametric  design  study  to  access  the  impact  of 
radome  shape,  source  location  and  dielectric  profiles  on  the 
radiation  pattern. 

4.  Perform  an  exhaustive  study  of  the  convergence  properties 
of  the  MM  solution. 

5.  To  improve  the  execution  time  some  of  the  Gaussian 
integration  loops  can  be  replaced  with  a  delta  function 
approximation  for  the  integral. 

6.  Symmetry  between  the  positive  and  negative  modes  can  also 
be  exploited  to  reduce  run  time. 

This  thesis  has  developed  flexible,  modular  code  for 
modeling  the  electric  field  radiation  pattern  for  a  source 
enclosed  in  a  dielectric  body  of  revolution.  Further  work 
should  concentrate  on  validation  and  improvements  to  reduce 
execution  time.  After  these  enhancements  are  incorporated,  the 
code  can  be  applied  to  the  radome  design  process. 


65 


PHERICAL  RADOME  (*),  AHALYT 


Figure  5.7  Effects  of  Radome  on  Far  Field  Radiation 
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APPENDIX  A.  SOURCE  CODE  AND  DATA  FILE 


A.  SOURCE  CODE  FOR  MAIN  PROGRAM 


C  PROGRAM  :  ldbormm.f 

C  DATE  :  23  January  1992 

C  REVISED  :  21  February  1992 

C  PROGRAMMERS  :  D.  JENN,  R.  FRANCIS 
C 

C  »  SPECIALIZED  FOR  ARBITRARY  CIRCULAR  APERTURE  EXCITATION  « 

C  »  EXCITATION  IS  SPECIFIED  IN  THE  SUBROUTINE  GENEX  (....)  « 
C 

c  BASED  ON  MAUTZ  AND  HARRINGTON'S  COMPUTER  PROG  FOR  BORS . 
c  all  or  part  of  the  surface  may  be  covered  with  a  surface 
c  impedance. 

C  SOURCE  IS  AT  BOR  COORDINATE  SYSTEM  ORIGIN. 

C 

c  imp=0  perfect  conductor 
c  iprint=0  print  pattern  points  to  unit  8 
c  iseg=0  print  the  generating  curve  points  to  unit  8 


C  BEGIN  MAIN  PROGRAM************ 

CHARACTER* 8  GNT , GNPHI , CGR, CGP 
CHARACTER* 14  TPTS , PPTS , PHIPTS , RHOPTS 
COMPLEX  EP,ET,Z( 100000) ,R(1600) ,B(1600) ,C(800) ,U,UC 
COMPLEXET1 , EP1 , ET2 , EP2 , EC, EX , ZLO (400) ,ZL(2400) , ETF , EPF 
COMPLEX  EXP1 , EXP2 , CONJG , CEXP, CMPLX , CT (3000) ,IMP,JX 
COMPLEX  CIRCTHETA, CIRCPHI 

DIMENSION  RH(400) ,ZH(400) ,XT(4) ,AT(4) ,IPS(800) 
DIMENSION  A ( 100) ,X(100) ,EXP(500) ,ANG(500) ,ECP(500) 
DIMENSION  XP(100) ,AP(100) ,XR(50) ,AR(50) 

DIMENSION  ECV (500) ,EXV(500) ,PHC(500) ,PHX(500) 

REAL  ETSCAT ( 500 ) ,EPSCAT(500) ,ETHF(500) ,EPHF(500) 
INTEGER  NT , NPHI , CNRHO , CNPHI , NP , SELECTION 
REAL  MODE, BASE, RS , ZP, RHB, ZHB 
DATA  PI, START, STOP/3. 1415926,0. ,90./ 

DATA  I PR I NT/ 0/ 

Rad=PI/180 . 

ECX=0. 

BK=2 . *PI 
U= ( 0 . , 1 . ) 
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U0=(0. ,0.) 

UC=-U/4./PI 

JK=( 0.0, 6. 283 185307) 

C  Select  the  subroutine  to  generate  BOR. 

WRITE (6, *) 'MAKE  BOR  SELECTION' 

WRITE (6,*)'  ' 

WRITE ( 6, *) 'Enter  number  from  menu  to  make  selection' 
WRITE (6 , *) '  ' 

WRITE ( 6 , * )  '  Selection  Number  BOR  Geometry' 

WRITE (6,*)'  1  OGIVE' 

WRITE (6,*)'  2  SPHERE' 

READ ( 5 , * ) SELECTION 

I F ( SELECTION . EQ . 1 ) THEN 

CALL  OGIVE (NP, ZH , RH , BASE , RS , ZP) 

ELSEIF (SELECTION . EQ . 2 ) THEN 

CALL  TESTS PHERE (NP , ZH , RH , BASE , RS , ZP) 

ENDIF 


C*********CALL  OGIVE  or  TESTSPHERE. 

IF  (NP.GT. 399) THEN 

WRITE (6,*) 'MAXIMUM  NUMBER  OF  POINTS (NP)  IS  399' 
GOTO  999 
ENDIF 

DT=90 . 0/ FLOAT (NP-1) 

WRITE (6,*) 'ENTER  THE  FILENAMES  gaus###' 

WRITE  (6,  *)  ' . ' 

WRITE ( 6, *) 'ENTER  THE  FILENAME  IN  T  (NT)' 

READ ( 5 , * ) GNT 

WRITE(6,*) 'ENTER  THE  FILENAME  IN  PHI  (NPHI) ' 

READ ( 5 , * ) GNPHI 

WRITE (6, *) 'ENTER  THE  FILENAME  IN  CNRHO' 

READ(5, *) CGR 

WRITE (6, *) 'ENTER  THE  FILENAME  IN  CNPHI' 

READ ( 5 , * ) CGP 

C  OPEN  THE  FILES  FOR  THE  gaus/gaus### 

TPTS=' gaus/ '//GNT 
PPTS= 'gaus/ '//GNPHI 
RHOPTS= ' gaus/ '//CGR 
PHIPTS='gaus/'//CG? 

OPEN  ( 1 ,  FILE=TPi‘S  ,  STATUS= '  OLD  '  ) 

OPEN ( 2 , FI LE=PPTS , STATUS= ' OLD ' ) 

OPEN ( 3 , FI LE=RHOPTS , STATUS^ ' OLD ' ) 

OPEN ( 4 , FILE=PHIPTS , STATUS= ' OLD' ) 

READ ( 1 , *) NT 

IF (NT. GT. 4) THEN 

WRITE ( 6, *) 'MAXIMUM  NUMBER  OF  POINTS (NT)  IS  4' 
GOTO  999 
ENDIF 

READ (2,*) NPHI 

IF (NPHI. GT. 2 00) THEN 
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WRITE ( 6 , * ) ' MAXIMUM  NUMBER  OF  POINTS(NPHI)  IS  200' 
GOTO  999 
ENDIF 

READ ( 3 , * ) CNRHO 

IF (CNRHO.GT. 50) THEN 

WRITE (6,*) 'MAXIMUM  NUMBER  OF  POINTS (CNRHO)  IS  50' 
GOTO  999 
ENDIF 

READ ( 4 , * ) CNPHI 

IF (CNPHI . GT. 100) THEN 

WRITE ( 6, *) 'MAXIMUM  NUMBER  OF  POINTS ( CNPHI )  IS  100' 
GOTO  999 
ENDIF 

C  LOAD  THE  WEIGHTS  AND  ABSCISSAS  IN  THE  VECTORS. 

DO  1  M=1 , NT 

READ (1 , * , END=1) XT (M) ,AT(M) 

1  CONTINUE 

DO  2  M=1 , NPHI 

READ ( 2  ,  *  , END=2 )  X  (M)  ,  A  (M) 

2  CONTINUE 

DO  3  M=l, CNRHO 

READ ( 3 , * , END=3 ) XR (M) , AR (M) 

3  CONTINUE 

DO  4  M=l, CNPHI 

READ (4, *,END=4)XP(M) ,AP(M) 

4  CONTINUE 
CLOSE (1) 

CLOSE (2) 

CLOSE (3) 

CLOSE (4) 

MP=NP-1 

MT=MP-1 

N=MT+MP 


C 


WRITE (6, *) 'ENTER 
READ ( 5, *) MODE 
WRITE (6, *) 'ENTER 
READ ( 5 , * ) P 
PHI=P*RAD 
WRITE (6,*) 'ENTER 
READ ( 5 , * ) SA 
SCAN=-SA*RAD 
WRITE (6, *) 'ENTER 
READ (5, *) IMP 
WRITE (6, *) 'ENTER 
READ ( 5 , * ) ARAD 


MODE  (FLOAT) ' 

PHI  (observation)  IN  DEGREES' 

THE  SCAN  ANGLE  IN  DEGREES' 

COMPLEX  IMPEDANCE' 

THE  ANTENNA  RADIUS  (wavelengths) ' 


ENTER  THE  NUMBER  OF  AZIMUTHAL  MODES 
(n=-MODES , . . . ,0, . . . ,+MODES) 

MODES=MODE 
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NBL0CK=2 *MODES+ 1 
MHI=M0DES+1 

OPEN ( 8 , FI LE= 7  out 1 dbor ' ) 

WRITE (8,8000)  2 . *BASE , NP, RS , ZP 
8000  FORMAT*//, 5X, '***  BOR  RADIATION  PATTERN  FOR  A  CIRCULAR ' 
*'  DISC  RADIATOR  USING  GENEX  ***' , 

*// , 2X, 'BOR  DIAMETER  (WAVELENGTHS) =  ' , F5 . 2 ,/, 2X , 

*  'NUMBER  OF  GENERATING  POINTS  (NP)=',I4, 

*// , 2X, 'SURFACE  RADIUS ', F5 . 2 , 2X, ' ZPRIME ' , F5 . 2 ) 

WRITE (8,30)  NT, NPHI 

30  FORMAT (/ , 12X, '  NT  NPHI ' ,/ , 10X, 15 , 2X, 15) 

IF(ISEG.EQ.O)  WRITE(8, 1300) 

1300  FORMAT (/ , 10X, 'INDEX', 8X, 'Z(I) ',10X, 'RHO(I) ',12X, 'SURF 

IMPED') 

DO  52  1=1, NP 

IF(ABS(ZH(I) ) .LT. .001)  ZH(I)=0. 

IF(ABS(RH(I) ) .LT. .001)  RH(I)=0. 

ZHB=ZH(I)/BK 
RHB=RH(I)/BK 
ZlO (i) =IMP 

ASSIGN  SURFACE  IMPEDANCE  AT  THIS  POINT.  THE  SURFACE 
IMPEDANCE  OF  SEGMENT  I  IS  ZLO(I) 

WRITE ( 8 , 8004 )  I , ZHB , RHB, ZLO ( I ) 

52  CONTINUE 

8004  FORMAT (11X,I4,4X,F8.3,8X, F8 . 3 , 6X, 2F8 . 2 ) 

C 

C  MODE  LOOP  TO  CALCULATE  THE  CURRENT  COEFFICIENTS.  POS  AND  NEG 
C  MODES  DONE  IN  THE  SAME  ITERATION  OF  THE  LOOP 


C*************ZLOAD,ZMAT, GENEX, DECOMP, SOLVE 

IF ( CABS ( IMP) . NE . 0 )  CALL  ZLOAD (NP, RH , ZH , ZLO , ZL) 

DO  400  M=1 , MHI 

NM=M-1 

CALL  ZMAT (NM, NM, NP, NPHI , NT , RH , ZH, X , A, XT, AT , Z ) 

IF (CABS (IMP) .NE.O)  CALL  ZTOT (MT, MP, ZL, Z ) 

CALL  GENEX ( NM , NP , NT , NPHI , CNRHO , CNPHI , XT , AT , X , A , 
*  XR,AR,XP,AP, SCAN, PHI, ARAD, RH,ZH,B) 

CALL  DECOMP (N, IPS, Z) 

CALL  SOLVE (N, IPS, Z,B,C) 


£★***********★**★*****★★**★★**★***★★** 

C  STORE  CURRENT  COEFFICIENTS  IN  ONE  LONG  COLUMN  VECTOR 
NTOPl=MODES-NM 
NTOP2=NBLOCK- (NTOP1+1) 

NS2=NT0P1*N 
NS1=NT0P2*N 
C  POSITIVE  MODE 
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DO  401  L=1,N 
401  CT (L+NS1) =C ( L) 

IF(NM.NE.O)  THEN 
C  NEGATIVE  MODE 
NMN=-NM 


C************* ******** *ZMAT,GENEX, DECOMP, SOLVE. 

CALL  ZMAT ( NMN , NMN , NP , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 

IF ( CABS ( IMP) . NE . 0 )  CALL  ZTOT (MT, MP, ZL, Z ) 

CALL  GENEX ( NMN , NP , NT , NPHI , CNRHO , CNPHI , XT , AT , X , A , 

*  XR,AR,XP,AP, SCAN, PHI, ARAD, RH,ZH,B) 

CALL  DECOMP (N, IPS, Z) 

CALL  SOLVE (N, IPS, Z,B,C) 
q* ********************************* 

DO  402  L=1,N 
402  CT (L+NS2 ) =C (L) 

ENDIF 

400  CONTINUE 
IT=NP 

DT=STOP/ FLOAT (NP-1 . ) 

BEGIN  PATTERN  LOOP  FROM  START  TO  STOP  IN  INCREMENTS  OF  DT 
(ALL  IN  DEG) 

DO  500  1=1, IT 
THETA=FLOAT ( 1-1 ) *DT+START 
THX=THETA*RAD 
PHR=PHR0 
RHB=RH(I)/BK 
ZHB=ZH(I)/BK 
I F ( THETA . GT .180.)  THEN 
PHR=PHR0+PI 
THX= (360. -THETA) *RAD 
ENDIF 

ET1= ( 0 . ,0. ) 

EP1= ( 0 . , 0. ) 

ET2= ( 0 . ,0. ) 

EP2= ( 0 . ,0. ) 

DO  300  M=1 , MHI 
NM=M-1 

EXP1=CEXP ( CMPLX ( 0 . , FLOAT (NM) *PHR) ) 

EXP2=CONJG (EXP1) 


c*****************plane 

CALL  PLANE (NM, NM, NP, NT, RH , ZH, XT, AT, THX, R) 

NTOPl=MODES-NM 

NTOP2=NBLOCK- (NTOP1+1) 

NS2=NT0P1*N 


73 


NS1=NT0P2*N 
DO  250  L=1,MT 

ET1=ET1+R(L) *CT ( L+NS1) *EXP1 
EP1=EP1+R ( L+N ) *CT(L+NS1) *EXP1 
IF(NM.EQ.O)  GO  TO  250 
ET1=ET1+R(L) *CT(L+NS2) *EXP2 
EP1=EP1-R(L+N) *CT(L+NS2) *EXP2 
250  CONTINUE 

DO  260  L=1,MP 

ET2=ET2+R(L+MT) *CT (L+NS1+MT) *EXP1 
EP2=EP2+R ( L+MT+N) *CT (L+NS1+MT) *EXP1 
IF(NM.EQ.O)  GO  TO  260 
ET2=ET2-R(L+MT) *CT (L+NS2+MT) *EXP2 
EP2=EP2+R (L+MT+N) *CT (L+NS2+MT) *EXP2 
260  CONTINUE 
300  CONTINUE 

C  DISK  CONTRIBUTION  IN  THE  FAR  FIELD  IS  ETF, EPF 
RRR=SQRT(RHB**2+ZHB**2) 

C*************CIRCTHETA, CIRCPHI 

ETF=CIRCTHETA ( CNPHI , XP , AP , CNRHO , XR , AR , ARAD , 

*  SCAN, PHI, RHB,ZHB) 

EPF=CIRCPHI ( CNPHI, XP,AP, CNRHO, XR, AR, ARAD, 

*  SCAN, PHI, RHB,ZHB) 

c  i  deleted  the  rrr*cexp ( jk*rrr)  factor  from  ETF  &  EPF. 

C* ****************************** ****** 

ETHF ( I ) =CABS ( ETF ) 

EPHF ( I ) =CABS (EPF ) 

C  TOTAL  E-THETA  AND  E-PHI  COMPONENTS 
ET=ET3  +ET2+ETF 
EP=EP1+EP2+EPF 
ETSCAT ( I ) =CABS ( ET-ETF) 

EPSCAT ( I ) =CABS ( EP-EPF) 

EC=ET 

EX=EP 

ECV ( I ) =CABS ( EC ) 

EXV ( I ) =CABS ( EX ) 

ECR=REAL(EC) 

ECI=AIMAG (EC) 

EXR=REAL(EX) 

EXI=AIMAG (EX) 

PHC ( I ) = AT AN 2 ( ECI , ECR+ 1 . e- 1 0 ) /RAD 
PHX ( I ) = AT AN 2 (EXI , EXR+1 . e-10) /RAD 
ANG ( I ) =THETA 

ECX=AMAX1 (ECX , ECV ( I ) ,EXV(I) ) 

500  CONTINUE 

WRITE (6,*)  'MAX  E  VALUE=/,ECX 
WRITE ( 8 , 103 )  P, ECX 

103  FORMAT (/, 1 OX, 'PHI  OF  RECEIVER  (DEG) = ' , F8 . 2 , / , 10X , 

*  'MAXIMUM  FIELD  VALUE  (V/M) =' , E15 . 5) 
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DO  600  1=1, IT 

ECP ( I ) = (ECV ( I ) /ECX) **2 

EXP ( I ) = (EXV ( I ) /ECX) **2 

ECP ( I ) =AMAX 1 ( ECP ( I ) , 1 . E-10) 

EXP ( I ) =AMAX1 ( EXP ( I ) , 1 . E-10) 

ECP (I ) =10 . *ALOG10 (ECP ( I ) ) 

EXP (I ) =10 . *ALOG10 (EXP (I ) ) 

600  CONTINUE 

IF ( IPRINT . EQ . 0)  THEN 
WRITE (8, 50 15) 

5015  FORMAT ( // , 7 X , 'ANGLE', 15X, 'CO-POL' , 25X, ' X-POL' , / , 7X , 

1' (DEG) ' ,4X,2(' (VOLTS) ' ,4X, ' (DEG) ' ,3X, ' (DB-REL) ' , 4X) ) 

DO  9000  L=1 , IT 

WRITE (8, 5016) ANG(L) ,ECV(L) ,PHC(L) ,ECP(L) ,EXV(L) ,PHX(L) 
1 , EXP (L) 

5016  FORMAT ( 5X ,F6.2,3X,2(F8.4,3X,F7.2,3X,F7.2,3X) ) 

9000  CONTINUE 

ENDIF 

OPEN  (2  ,  f  ile='bang.in' ) 

OPEN (3 , f ile='bcpole .m' ) 

OPEN (4 , f ile='bxpole.m' ) 

OPEN (7 , f ile='etscat .m' ) 

OPEN (8 , f ile='epscat .m' ) 

OPEN (9, file='etf .m' ) 

OPEN ( 10 , f ile= ' epf . m' ) 

DO  9097  1=1, IT 
WRITE (2, 5019)  ANG(I) 

WRITE (3 , 5019)  ECV (I) 

WRITE (4, 5019)  EXV (I) 

C  change  ecv,ecv  to  ecp  &  exp  to  get  normalized  db  values. 
WRITE (7, 5019)  ETSCAT(I) 

WRITE (8, 5019)  EPSCAT(I) 

WRITE (9, 5019)  ETHF(I) 

9097  WRITE (10, 5019)  EPHF(I) 

5019  FORMAT (F8 . 3 ) 

CLOSE (2) 

CLOSE (3) 

CLOSE (4) 

CLOSE (7) 

CLOSE (8) 

CLOSE (9) 

CLOSE ( 10 ) 

999  CONTINUE 
STOP 
END 

C* *********** **END  OF  MAIN  PROGRAM. 
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c 

C**************SUBROUTINE  ZMAT. 

C  REFERENCE:  AN  IMPROVED  E-FIELD  SOLUTION  FOR  CONDUCTING  BOR 
C  J.R.MAUTZ  AND  R.F.  HARRINGTON 

C  TECHNICAL  REPORT  TR-80-1 

C  ROME  AIR  DEVELOPMENT  CENTER 

C  CONTRACT  NO.  F-30602-79-C-0011 

SUBROUTINE  ZMAT (Ml , M2 , NP , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 

COMPUTE  THE  MM  IMPEDANCE  MATRIX  ELEMENTS.  THIS  IS  FROM 
MAUTZ  AND  HARRINGTON  (NO  CHANGES) . 


COMPLEXZ (100000) , U1 , U2 , U3 , U4 , U5 , U6 , U7 , U8 , U9 , UA, UB, G4A ( 4 ) ,G5A 
(4) 

COMPLEXG6A ( 4 )  ,G4B(4)  ,G5B(4)  ,G6B(4)  , H4A, H5A, H6A , H4  B , H5B 
COMPLEX  H6B, UC , UD, GA ( 100 ) ,GB(100) 

DIMENSION  RH (400) ,ZH(400) ,X(100) ,A(100) ,XT(4) ,AT(4) 
DIMENSION  RS (400) ,ZS(400) ,D(400) ,DR(400) ,DZ(400) 
DIMENSION  DM (400) ,  C2 (100) ,C3(100) ,R2(4) ,Z2(4) 

DIMENSION  C4 ( 100) ;C5(100) ,C6(100) , Z7 (4 ) , R7 ( 4) , Z8 (4 ) , R8 ( 4 ) 
CT=2  . 

CP=.  1 

DO  10  1=2, NP 
12=1-1 

RS ( 1 2 )  = . 5  * ( RH ( I ) +RH (12)) 

ZS (12) =. 5* (ZH (I) +ZH ( 12 ) ) 

Dl=. 5* (RH ( I ) -RH ( 12 ) ) 

D2= . 5* ( ZH ( I ) -ZH ( 12 ) ) 

D ( 12 ) =SQRT (D1*D1+D2*D2) 

DR ( 12 ) =D1 
DZ (12) =D2 

I F ( RS (12). EQ . 0 . )  RS ( 12 ) =1 . 

DM ( 12 ) =D ( 12 ) /RS ( 12 ) 

10  CONTINUE 
M3=M2-M1+1 
M4=M1-1 
PI2=1. 570796 
DO  11  K=1 , NPHI 
PH=PI2* (X (K) +: 

C2 (K) =PH*PH 
SN=SIN ( . 5*PH 
C3 (K) =4 . *SN*SN 
A1=FI2*A (K) 

D4=. 5*A1*C3 (K) 

D5=A1*C0S ( PH) 

D6=A1*SIN ( PH) 

M5=K 

DO  29  M=1 , M3 
PHM= (M4+M) *  PH 
A2=COS ( PHM) 
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C4 (M5) =D4*A2 
C5(M5)=D5*A2 
C6 (M5) =D6*SIN (PHM) 
M5=M5+NPHI 
29  CONTINUE 

11  CONTINUE 
MP=NP-1 
MT=MP-1 
N=MT+MP 
N2N=MT*N 
N2=N*N 
U1=(0. , .5) 

U2= ( 0 . ,2.) 

JN=-1-N 

DO  15  JQ=1 , MP 
KQ=2 

IF(JQ.EQ.l)  KQ=1 
IF(JQ.EQ.MP)  KQ=3 
R1=RS ( JQ) 

Z1=ZS (JQ) 

D1=D ( JQ) 

D2=DR ( JQ) 

D3=DZ (JQ) 

D4=D2/R1 
D5=DM ( JQ) 

SV=D2/D1 

CV=D3/D1 

T6=CT*D1 

T62=T6+D1 

T62=T62*T62 

R6=CP*R1 

R62=R6*R6 

DO  12  L=1 , NT 

R2 (L)=R1+D2*XT(L) 

Z2 (L)=Z1+D3*XT(L) 

12  CONTINUE 
U3=D2*U1 
U4=D3*U1 

DO  16  IP=1,MP 
R3=RS (IP) 

Z3=ZS (IP) 

R4=R1-R3 
Z4=Z 1-Z  3 
FM=R4*SV+Z4*CV 
PHM=ABS(FM) 

PH=ABS (R4*CV-Z4*SV) 
D6=PH 

IF(PHM.LE.Dl)  GO  TO  26 
D6=PHM-D1 

D6=SQRT(D6*D6+PH*PH) 

26  IF ( IP. EQ . JQ)  GO  TO  27 


77 


KP=1 

IF (T6 .GT. D6)  KP=2 
IF(R6.GT.D6)  KP=3 
GO  TO  28 

27  KP=4 

28  GO  TO  (41,42,41,42) , KP 
42  DO  40  L=1 , NT 

D7=R2 (L) -R3 
D8=Z2 (L) -Z3 
Z7 (L) =D7*D7+D8*D8 
R7 (L)=R3*R2 (L) 

Z8(L)  =  .25*Z7(L) 

R8(L)=.25*R7(L) 

40  CONTINUE 

Z4=R4*R4+Z4*Z4 
R4=R3*R1 
R5= . 5*R3*SV 
DO  33  K=1 , NPHI 
A1=C3 (K) 

RR=Z4+R4*A1 
UA=0 . 

UB=0 . 

IF (RR. LT.T62)  GO  TO  34 
DO  35  L=1 , NT 
R=SQRT(Z7 (L) +R7 (L) *A1) 

SN=-SIN (R) 

CS=COS (R) 

UC=AT ( L) /R*CMPLX ( CS , SN) 
UA=UA+UC 
UB=XT ( L) *UC+UB 
35  CONTINUE 
GO  TO  36 
34  DO  37  L=1 , NT 

R=SQRT ( Z8 ( L) +R8 ( L) *A1) 
SN=-SIN(R) 

CS=COS (R) 

UC=AT ( L) /R*SN*CMPLX ( -SN , CS ) 

UA=UA+UC 

UB=XT(L) *UC+UB 

37  CONTINUE 
A2=FM+R5*A1 
D9=RR-A2  *A2 
R=ABS ( A2 ) 

D7=R-D1 

D8=R+D1 

D6=SQRT(D8*D8+D9) 

R=SQRT(D7*D7+D9) 

IF ( D7 . GE . 0 . )  GO  TO  38 
Al=ALOG ( (D8+D6) * (-D7+R) /D9) /D1 
GO  TO  39 

38  Al=ALOG( (D8+D6)/ (D7+R) )/Dl 
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39  UA=A1+UA 

UB=A2  * ( 4 . / ( D6+R) -A1 ) /Dl+UB 
36  GA(K)=UA 
GB (K) =UB 
33  CONTINUE 
K1=0 

DO  45  M=1 , M3 
H4A=0 . 

H5A=0 . 

H6A=0 . 

H4B=0. 

H5B=0. 

H6B=0 . 

DO  46  K=1,NPHI 
K1=K1+1 
D6=C4 (Kl) 

D7=C5(K1) 

D8=C6 (Kl) 

UA=GA(K) 

UB=GB(K) 

H4A=D6*UA+H4A 
H5A=D7*UA+H5A 
H6A=D8*UA+H6A 
H4B=D6*UB+H4B 
H5B=D7*UB+H5B 
H6B=D8*UB+H6B 
46  CONTINUE 
G4A(M)=H4A 
G5A(M)«H5A 
G6A (M) =H6A 
G4B (M) =H4B 
G5B (M) =H5B 
G6B (M) =H6B 
45  CONTINUE 

IF (KP. NE . 4 )  GO  TO  47 
A2=D1/ (PI2*R1) 

D6=2./D1 
D8=0 . 

DO  63  K=1 , NPHI 
A1=R4  *C2 (K) 

R=R4  *C3 (K) 

IF (R. LT.T62)  GO  TO  64 
D7=0. 

DO  65  L=1 , NT 

D7=D7+AT(L)/SQRT(Z7 (L)+A1) 

65  CONTINUE 
GO  TO  66 

64  A1=A2/(X(K)+1. ) 

D7=D6*ALOG (Al+SQRT ( 1 . +A1*A1) ) 

66  D8=D8+A(K) *D7 
63  CONTINUE 
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Al= . 5*A2 
A2=l ./Al 

D8=-PI2*D8+2 . /Rl* (BLOG ( A2 ) +A2*BLOG (Al) ) 

DO  67  M=1 , M3 
G5A (M) =D8+G5A (M) 

67  CONTINUE 
GO  TO  47 
41  DO  25  M=1,M3 
G4A (M) =0 . 

G5A (M) =0 . 

G6A (M) =0 . 

G4B (M) =0 . 

G5B (M) =0 . 

G6B (M) =0 . 

25  CONTINUE 

DO  13  L=1 , NT 
A1=R2 (L) 

R4=A1-R3 

Z4=Z2 (L) -Z3 

Z4=R4*R4+Z4*Z4 

R4=R3*A1 

DO  17  K=1 , NPHI 

R=SQRT (Z4+R4*C3 (K)  ) 

SN=-SIN(R) 

CS=COS(R) 

GA (K) =CMPLX (CS  f  SN) /R 
17  CONTINUE 
D6=0 . 

IF (R62 . LE . Z4 )  GO  TO  51 

DO  62  K=1 , NPHI 

D6=D6+A (K) /SQRT ( Z4+R4  *C2 (K) ) 

62  CONTINUE 

Z4=3 . 14 1593/SQRT ( Z4/R4 ) 

D6=-PI2*D6+ALOG ( Z4+SQRT (l.+Z4*Z4) )/SQRT(R4) 
51  A1=AT (L) 

A2=XT(L) *A1 
K1=0 

DO  30  M=1 , M3 
U5=0 . 

U6=0 . 

U7=0. 

DO  32  K=1 , NPHI 
UA=GA (K) 

K1=K14 1 

U5=C4 (Kl) *UA+U5 
U6=C5 (Kl) *UA+U6 
U7=C6 (Kl) *UA+U7 
32  CONTINUE 
U6=D6+U6 

G4A (M) =A1*U5+G4 A (M) 

G5A (M) =A1*U6+G5A (M) 
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G6A (M) =A1*U7+G6A (M) 
G4B (M) =A2*U5+G4B (M) 
G5B(M)=A2*U6+G5B(M) 
G6B(M)=A2*U7+G6B(M) 

30  CONTINUE 

13  CONTINUE 

47  A1=DR ( IP) 

UA=A1*U3 
UB=DZ ( IP) *U4 
A2=D ( IP) 

D6=-A2*D2 

D7=D1*A1 

D8=D1*A2 

JM=JN 

DO  31  M=1 , M3 
FM=M4+M 
A1=FM*DM (IP) 

H5A=G5A (M) 

H5B=G5B (M) 

H4A=G4A (M) +H5A 
H4B=G4B (M) +H5B 
H6A=G6A (M) 

H6B=G6B(M) 

U7=UA*H5A+UB*H4A 

U8=UA*H5B+UB*H4 B 

U5=U7-U8 

U6=U7+U8 

U7—-Ui*H4A 

U8=D6*H6A 

U9=D6*H6B-A1*H4A 

UC=D7* (H6A+D4  *H6B) 

UD=FM*D5*H4A 

K1=IP+JM 

K2=K1+1 

K3=K1+N 

K4=K2+N 

K5=K2+MT 

K6=K4+MT 

K7=K3+N2N 

K8=K4+N2N 

GO  TO  (18, 20, 19) ,KQ 

18  Z(K6)=U8+U9 
IF(IP.EQ.l)  GO  TO  21 
Z(K3)=Z(K3) +U6-U7 
Z(K7)=Z(K7) +UC-UD 
IF(IP.EQ.MP)  GO  TO  22 

21  Z  (K4 )  =U6-(-U7 
Z (K8 ) =UC+UD 
GO  TO  22 

19  Z (K5) =Z (K5) +U8-U9 
IF(IP.EQ.l)  GO  TO  23 
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non 


Z  (K1)=Z  (Kl)  +U5+L  7 
Z ( K7 ) =Z ( K7 ) +UC-UD 
IF (IP. EQ.MP)  GO  TO  22 

23  Z (K2 ) =Z (K2 ) +U5-U7 
Z (K8)=UC+UD 

GO  TO  22 

20  Z (K5) =Z (K5) +U8-U9 
Z(K6)=U8+U9 
IF(IP.EQ.l)  GO  TO  24 
Z (K1)=Z (Kl) +U5+U7 
Z (K3 ) =Z (K3 ) +U6-U7 
Z (K7 ) =Z (K7 ) +UC-UD 
IF(IP. EQ.MP)  GO  TO  22 

24  Z (K2) =Z (K2) +U5-U7 
Z (K4)=U6+U7 
Z(K8)=UC+UD 

22  Z (K8+MT) =U2* (D8* (H5A+D4*H5B) -A1*UD) 
JM=JM+N2 

31  CONTINUE 

16  CONTINUE 
JN=JN+N 

15  CONTINUE 
RETURN 
END 


C************* * *SUBROUTINE  SOLVE . 

C  REFERENCE  :  TECHNICAL  REPORT  TR-80-1 
SUBROUTINE  SOLVE (N , IPS , UL, B, X) 

SEE  MAUTZ  &  HARRINGTON  FOR  DETAILS 

COMPLEX  UL(  100000)  ,  B  (  800)  ,  X  (800)  . 
DIMENSION  IPS (800) 

NP=N+1 
IP=IPS(1) 

X  ( 1 )  =B  (IP) 

DO  2  1=2, N 
IP=IPS(I) 

IPB=IP 
IM1=I-1 
SUM= ( 0 . , 0 . ) 

DO  1  J=1 , IM1 
SUM=SUM+UL (IP) *X(J) 

1  IP=IP+N 

2  X(I)=B(IPB) -SUM 
K2=N* (N-l ) 

IP=IPS (N) +K2 
X (N) =X (N) /UL (IP) 

DO  4  IBACK=2 , N 
I=NP-IBACK 
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non 


K2=K2-N 
IPI=IPS (I) +K2 
IP1=I+1 
SUM= ( 0 . ,0. ) 

IP=IPI 

DO  3  J=IP1 , N 
IP=IP+N 

3  SUM=SUM+UL ( IP) *X(J) 

4  X(I)=(X(I) -SUM)/UL(IPI) 
RETURN 

END 


C*  *  *  *  *  *SUBROUTINE  DECOMP 
C  REFERENCE:  TECHNICAL  REPORT  TR-80-1 
SUBROUTINE  DECOMP (N, IPS , UL) 

SEE  MAUTZ  &  HARRINGTON  FOR  DETAILS 


COMPLEX  UL(IOOOOO) , PIVOT, EM 
DIMENSION  SCL(800) , IPS (800) 

DO  5  1=1, N 
IPS(I)=I 
RN=0 . 

J1=I 

DO  2  J=1 , N 

ULM=ABS(R£AL(UL(J1) ) ) +ABS (AIMAG (UL(J1) ) ) 

J1=J1+N 

IF (RN-ULM)  1,2,2 

1  RN=ULM 

2  CONTINUE 
SCL ( I ) =1 . /RN 

5  CONTINUE 
NM1=N-1 
K2=0 

DO  17  K=1 , NM1 
BIG=0 . 

DO  11  I=K,N 
IP=IPS(I) 

IPK=IP+K2 

SIZE= (ABS (REAL(UL(IPK) )) +ABS (AIMAG (UL ( IPK) ) ) ) *SCL(IP) 
IF ( SIZE-BIG)  11,11,10 

10  BIG=SIZE 
IPV=I 

11  CONTINUE 

IF (IPV-K)  14,15,14 

14  J=IPS (K) 

IPS (K) =IPS ( IPV) 

IPS(IPV)=J 

15  KPP=IPS (K) +K2 
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o  o  o  o 


PIVOT=UL(KPP) 

KP1=K+1 
DO  16  I=KP1 , N 
KP=KPP 
IP=IPS (I) +K2 
EM=-UL( IP) /PIVOT 
18  UL(IP) =-EM 
DO  16  J=KP1,N 
IP— IP+N 
KP=KP+N 

UL(IP) =UL(IP) +EM*UL(KP) 

16  CONTINUE 
K2=K2+N 

17  CONTINUE 
RETURN 
END 


C************FUNCTION  BLOG 

C  REFERENCE:  TECHNICAL  REPORT  TR-80-1  ; (page  56) 
FUNCTION  BLOG (X) 

IF (X. GT. . 1)  GO  TO  1 
X2=X*X 

BLOG= ( ( .075*X2-. 1666667) *X2+1. ) *X 
RETURN 

1  BLOG=ALOG (X+SQRT ( 1 . +X*X) ) 

RETURN 

END 


C* ***************  *  *SUBROUTINE  PLANE 

C  REFERENCE:  TECHNICAL  REPORT  TR-80-1  ; (pages  57-64) 
SUBROUTINE  PLANE (Ml , M2 , NP , NT , RH , ZH , XT , AT , THR , R) 

PLANE  WAVE  EXCITATION  VECTOR.  FROM  MAUTZ  AND  HARRINGTON. 
MODIFIED  TO  DO  ONLY  ONE  ANGLE  AND  FREQUENCY  PER  CALL. 

COMPLEX  R ( 1600) , U , U1 , UA , UB, FA ( 5) ,FB(5) , F2A, F2B, FI A, FIB 
COMPLEX  U2 , U3 ,  U4 , U5 

DIMENSION  RH (400) ,ZH(400) , XT ( 4 ) , AT (4 ) , R2 (4 ) 

DIMENSION  Z2 ( 4 ) , BJ ( 2000) 

MP=NP-1 
MT=MP-1 
N=MT+MP 
N2=2*N 
CC=COS (THR) 

SS=SIN (THR) 

U=(0. , 1. ) 

U1=3.141593*U**M1 
M3=M1+1 
M4=M2+3 
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IF(Ml.EQ.O)  M3=2 

M5=Ml+2 

M6=M2+2 

DO  12  IP=1;MP 

K2=IP 

I»IP+1 

DR= . 5  * (RH(I) -RH (IP) ) 

DZ=. 5* (ZH (I) -ZH (IP) ) 

D1=SQRT ( DR*DR+DZ  *DZ ) 

Rl= .25* (RH ( I ) +RH (IP) ) 

IF (R1 . EQ . 0 . )  Rl=l. 

Zl= . 5* ( ZH ( I ) +ZH (IP) ) 

DR=.5*DR 

D2=DR/R1 

DO  13  L=1 , NT 

R2 (L)=R1+DR*XT(L) 

Z2 (L)=Z1+DZ*XT(L) 

13  CONTINUE 
D3=DR*CC 
D4=-DZ*SS 
D5=D1*CC 
DO  23  M=M3 , M4 
FA (M) =0 . 

FB (M) =0 . 

23  CONTINUE 

DO  15  L=1 , NT 
X=SS*R2 (L) 

IF(X.GT. .5E-7)  GO  TO  19 
DO  20  M=M3 , M4 
BJ (M) =0 . 

20  CONTINUE 
BJ  ( 2 )  =1 . 

S=1 . 

GO  TO  18 

19  M=2.8*X+14.-2./X 

I F ( X . LT . . 5 )  M=ll . 8+ALOG10 (X) 
IF (M. LT.M4)  M=M4 
BJ  (M)  =0 . 

JM=M- 1 
BJ ( JM) =1 . 

DO  16  J=4 , M 
J2=JM 
JM=JM-1 
J1=JM-1 

BJ ( JM) =J1/X*BJ ( J2 ) -BJ ( JM+2 ) 

16  CONTINUE 
S=0. 

IF (M . LE . 4 )  GO  TO  24 
DO  17  J=4,M,2 
S=S+BJ (J) 

17  CONTINUE 
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24  S=BJ (2 ) +2 . *S 
18  ARG=Z2(L)*CC 

UA=AT (L) /S*CMPLX (COS (ARG) , SIN (ARG) ) 
UB=XT (L) *UA 
DO  25  M=M3,M4 
FA  (M)  =BJ  (M)  *UA+FA  (M) 

FB (M) =BJ (M) *UB+FB (M) 

25  CONTINUE 
15  CONTINUE 

IF(Ml.NE.O)  GO  TO  26 
FA(1) =FA(3) 

FB ( 1 ) =FB ( 3 ) 

26  UA=U1 

DO  27  M=M5,M6 

M7=M-1 

M8=M+1 

F2A=UA* ( FA (M8 ) +FA (M7 ) ) 

F2B=UA* (FB(M8)+FB(M7) ) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(M7) ) 
F1B=UB*(FB(M8)-FB(M7) ) 

U4=D4*UA 

U2=D3*F1A+U4*FA(M) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R(K2+MT) =-D5* (F2A+D2*F2B) 

R (K5+MT) =D1* (F1A+D2*F1B) 

IF(IP.EQ.l)  GO  TO  21 
R(K1) =R(K1) +U2-U3 
R (K4 ) =R (K4 ) +U4-U5 
IF(IP.EQ.MP)  GO  TO  22 

21  R(K2)=U2+U3 
R(K5)=U4+U5 

22  K2=K2+N2 
UA=UB 

27  CONTINUE 
12  CONTINUE 

RETURN 

END 


C*************SUBR0UTINE  ZLOAD. 

C  REFERENCE:  COMPUTATION  OF  RADIATION  AND  SCATTERING  FOR 
C  LOADED  BODIES  OF  REVOLUTION 

C  HARRINGTON  AND  MAUTZ 

C  REPORT  AFCRL-70-0046 

C  AIRFORCE  CAMBRIDGE  RESEARCH  LABS 

C  CONTRACT  NO.  F-19628-68-C-0180 


86 


SUBROUTINE  ZLOAD(NP,RH, ZH, ZO, Z) 


C 

C  COMPUTES  IMPEDANCE  MATRIX  ELEMENTS  FOR  LOADED  BODIES  OF  REV 
C  Z0(I)  IS  THE  SURF  IMPEDANCE  OF  THE  ITH  SEGMENT  (NP-1 
C  SEGMENTS)  Z(.)  ARE  THE  IMPEDANCE  MATRIX  TERMS  (TRIDIAGONAL 
C  FOR  T-T  SUBMATRIX;  DIAGONAL  FOR  P-P  SUBMATRIX) .  STORED  IN 
C  COL  VECTOR. 

C 


COMPLEXC1 , C2 , ZO (400) ,Z(2400) , XI , X2 , X3 , Y1 , Y2 , Y3 , FN ( 4 00) 
COMPLEX  U1 , U2 , U3 , XI , YI 

DIMENSION  RH (400) fZH(400) ,RS(400) ,D(400) ,SV(400) 

PI=3. 14159 

MT=NP-2 

MP=NP-1 

N=MT+MP 

DO  10  IP=2 , NP 

II=IP-1 

DR=RH ( IP) -RH(II) 

DZ=ZH (IP) -ZH ( II ) 

D ( II ) =SQRT (DR*DR+DZ  *DZ ) 

SV ( II ) =DR/D ( II ) 

RS ( II ) = . 5* (RH (IP) +RH ( II ) ) 

DS=D ( II ) *SV ( II ) /2 . 

Q1=RS ( II ) +DS 
Q2=RS(II) -DS 
FN ( II ) =1 . 

IF ( (ABS (Q2 ) .GT.l.E-6) .AND. (ABS(Ql) .GT.l.E-6) ) 

*  FN ( II ) =ALOG (Q1/Q2 ) 

10  CONTINUE 
L0=MT*3-2 
DO  20  1=1, MP 
C1=PI*Z0 ( I ) 

IF(I.EQ.MP)  GO  TO  80 
KI=2 

IF(I.EQ.l)  KI=1 
IF(I.EQ.MT)  KI=3 
11=1+1 
C2=PI*Z0 (II) 

A=SV ( I ) 

IF (ABS (A) .LT.l.E-6)  GO  TO  41 
Xl=Cl*FN(I)/2./A 

X2=C1*2 . / A* ( 1 . -RS ( I ) *FN (I)/D(I) /A) 

X3=-X2*RS(I)/D(I)/A 
GO  TO  42 

41  CONTINUE 
Xl=Cl/2./RS(I)*D(I) 

X2= ( 0 . ,0. ) 

X3=Cl*D(I)/6./RS(I) 

42  CONTINUE 
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A=SV ( II ) 

IF(ABS(A) .LT.l.E-6)  GO  TO  45 
Yl=C2*FN(II)/2./A 

Y2=C2*2 ./A* ( 1 . -RS ( II ) *FN ( II ) /D ( II ) /A) 

Y3=-Y2*RS (II)/D(II)/A 
GO  TO  40 
45  CONTINUE 

Yl=C2/2 ./RS (II)  *D (II) 

Y2= ( 0 . ,0.) 

Y3=C2*D(II) /6 ./RS (II) 

40  CONTINUE 
C 

C  DEFINE  TRIDIAGONAL  ELEMENTS  FOR  T-T  SUBMATRIX  (STORED  IN 
C  COLS)  (Ul-  DIAG ;  U2-  LOWER;  U3-  UPPER) 

C 

XI=X1+X2+X3 

YI=Y1-Y2+Y3 

IF(KI.EQ.l)  XI=C1/SV(I) 

C  IF (KI . EQ . 3 )  YI=C2/SV(II) 

U1=XI+YI 

U2=X1-X3 

U3=Y1-Y3 

L=2  + (1-2) *3 

IF(KI.EQ.l)  L=0 

L1=L+1 

L2=L+2 

L3=L+3 

go  to  (50,60,70) , ki 
50  Z ( LI ) =U1 
Z (L2) =U2 
GO  TO  8  0 
60  Z(L1)=U3 
Z (L2)=U1 
Z (L3) =U2 
GO  TO  80 
70  Z(L1)=U3 
Z(L2)=U1 

80  Z(L0+I)=2.*C1*D(I)/RS(I) 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  ZTOT (MT , MP , ZL, Z ) 

C 

C  ADDS  THE  SURF  IMPEDANCE  TERMS  TO  THE  TRIDIAGONAL  ELEMENTS  OF 
C  THE  BOR  IMPEDANCE  MATRIX  Z. 

C 

COMPLEX  ZL(2400) ,Z(100000) 

N=MT+MP 
M0=MT*3-2 
DO  100  1=1, MP 
L0=MT*N+ ( I — 1 ) *N+MT 
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IF(I.EQ.MP)  GO  TO  80 
KI=2 

IF(I.EQ.l)  KI=1 
IF(I.EQ.MT)  KI=3 
L2= ( I — 1 ) *N+I 
L1=L2 -1 
L3=L2+1 
M=2+3* (1-2) 

IF(KI.EQ.l)  M=0 

M1=M+1 

M2=M+2 

M3=M+3 

go  to  (50 , 60 , 70) , ki 
50  Z (L2)=Z (L2)+ZL(M1) 
Z(L3)=Z(L3)+ZL(M2) 

GO  TO  80 

60  Z(L1)=Z(L1)+ZL(M1) 
Z(L2)=Z ( L2 ) +ZL (M2 ) 

Z ( L3 ) =Z ( L3 ) +ZL (M3 ) 

GO  TO  80 

70  Z(L1)=Z(L1)+ZL(M1) 
Z(L2)=Z(L2)+ZL(M2) 

80  Z (LO+I ) =Z ( LO+I ) +ZL (MO+I ) 
100  CONTINUE 
RETURN 
END 


C*********SUBROUTINE  OGIVE 
C  SUBROUTINE  :  OGIVE 
C  DATE  :  4  SEPTEMBER  1991 

C  PROGRAMMER  :  R.M.  FRANCIS 
C  REVISED  :  26  JANUARY  1992 

C  COMMENTS  : 

C  THIS  SUBROUTINE  WILL  GENERATE  DATA  FOR  A  BODY  OF  REVOLUTION 
C  (BOR)  IN  THE  FORM  OF  AN  OGIVE. 

C  DIMENSIONS  ARE  NORMALIZED  TO  WAVELENGTH. 

C  ZH  =  Z  CO-ORDINATE  *  2*PI 
C  RH  =  RADIUS  *2*PI 

SUBROUTINE  OGIVE (NP , ZH , RH , BASE) 

C  NP  =  NUMBER  OF  POINTS  ON  THE  OGIVE  SURFACE,  MAXIMUM  =  1000 
C  ZP  =  Z PRIME ,  THE  POSITION  ON  Z  WHERE  THE  RADIUS  OF 
C  CURVATURE  STARTS 
C  BASE  =  BASE  RADIUS 

C  RS  =  RADIUS  OF  CURVATURE  IN  THE  RZ , Z  PLANE  OF  THE  OGIVE 
INTEGER  I , NP 

REAL  RH(400) ,ZH(400) , ZP, BASE , RS , ZCOORD, RADIUS , AL 
PI=3. 1415926 
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C  INPUT  THE  VARIABLES  FOR  THE  OGIVE, ZP, B, RS , NP 

WRITE (6,*)  ' ENTER  SURFACE  CURVATURE  (wavelengths)' 

READ ( 5 , * )  RS 

WRITE  ( 6 ,  * )  # ENTER  ZPRIME ,  WHERE  CURVATURE  STARTS  (wavelengths )  7 
READ ( 5 , * )  ZP 

WRITE (6 , *)  'ENTER  BASE  RADIUS  (wavelengths)' 

READ (5,*)  BASE 

C  PERFORM  CALCULATIONS 

AL=SQRT (2 . *BASE*RS-BASE**2 ) 

ZMAX=  AL  +  ZP 
ANG=ASIN (AL/RS) 

L=ANINT( (RS*ANG+ZP) *5) 

NP=2*L+1 

WRITE (6,*) 'NUMBER  OF  POINTS  IS:  ',NP 
DZ=  ZMAX/ FLOAT (NP-1) 

DO  10  1=1, NP 

ZCOORD=  (1-1) *DZ 
ZH ( I ) =2*PI*ZCOORD 

IF  (ZCOORD.LE.ZP)  THEN 
RADIUS=BASE 

ELSE 

RADIUS=SQRT (RS**2- ( ZCOORD-ZP) **2) + (BASE-RS) 
ENDIF 

RH (I) =2*PI*RADIUS 
10  CONTINUE 
RETURN 
END 


C*********SUBR0UTINE  GENEX. 

C  SUBROUTINE  :  GENEX 
C  DATE  :  14  JANUARY  1992 

C  REVISED  :  17  February  1992 

C  PROGRAMMER  :  R.M.  FRANCIS 

SUBROUTINE  GENEX ( NM , NP , NT , NPHI , CNRHO , CNPHI , XT , AT , X , A , 

*  XR,AR,XP,AP, SCAN, PHI, ARAD, RH,ZH,R) 

C  SOURCE  ON  Z  AXIS  AT  Z=0. 

C  THE  EXCITATION  VECTOR  IS  COMPUTED  FOR  THE  GIVEN  R, THETA  AND 
C  PHI  COMPONENTS  SPECIFIED  IN  FUNCTION  SUBROUTINES 
C  CIRCTHETA , CIRCRHO , CIRCPHI . 

COMPLEX  R (1600) , CEXP, PSI , ST1 , ST2 , ST3 , ST4 , SP, SI , S2 , S3 , S4 , S5 
COMPLEX  CIRCTHETA, CIRCRHO, CIRCPHI 

DIMENSION  RH (400) ,ZH(400) ,XT(4) ,AT(4) ,X(100) ,A(100) 
REAL  XP(100) , AP(100) ,XR(50) ,AR(50) , SCAN, PHI , ARAD, RHB, ZHB 
INTEGER  NM , NP , NT , NPHI , CNRHO , CNPHI , I , N , MP , MT 
PI=3. 141592654 


BK=2 . *PI 
MP=NP-1 
MT=MP-1 
N=MT+MP 

C  LIMITS  ON  PHI  INTEGRATION 
Pl=(2.*PI-0. )/2. 

P2=(2.*PI+0. )/2. 

DO  30  IP=1,MP 

C  QUANTITIES  FOR  THE  FIRST  SEGMENT  (POSITIVE  SLOPE) 

I=IP+1 

II=IP 

DR=RH ( I ) -RH(II) 

DZ=ZH(I) -ZH(II) 

D1=SQRT ( DR*DR+DZ *  DZ ) 

R1=.5*(RH(I)+RH(II) ) 

Zl=. 5* ( ZH (I) +ZH (II ) ) 

SVP1=DR/D1 

CVP1=DZ/D1 

V1=ATAN2 (DR, DZ+1 . E-5) 

C  QUANTITIES  FOR  THE  SECOND  SEGMENT  (NEGATIVE  SLOPE) 

C  (SKIP  IF  LAST  SEGMENT) 

I=IP+2 

II=IP+1 

DR=RH ( I ) -RH (II) 

DZ=ZH ( I ) -ZH ( II ) 

D2=SQRT ( DR*DR+DZ*DZ ) 

R2= . 5* (RH ( I ) +RH ( II ) ) 

Z2= . 5* (ZH (I) +ZH (II) ) 

SVP2=DR/D2 

CVP2=DZ/D2 

V2=ATAN2 (DR, DZ+l.E-5) 

C  BEGIN  PHI  INTEGRATION:  R  TERMS  IN  SI  AND  S2 ;  THETA  TERM  IN 
C  S3  AND  S4 ;  PHI  TERM  IN  S5. 

ST1= ( 0 . ,0. ) 

ST2= ( 0 . ,0. ) 

ST3=(0. ,0.) 

ST4=(0. ,0. ) 

SP= ( 0 . , 0 . ) 

DO  20  J=1 , NPHI 
PH=P1*X ( J) +P2 

PSI=CEXP(CMPLX(0. ,-MODE*PH) ) 

S1=(0. ,0. ) 

S2= ( 0 • ,0. ) 

S3= ( 0 • , 0. ) 

S4= ( 0 . ,0. ) 

S5=(0. ,0. ) 

IF(IP.LE.MT)  THEN 

C  t-CURRENT  INTEGRATION  FOR  THE  POSITIVE  SLOPE 
DO  13  L=1 , NT 
TP=D1*XT(L) /2. 

RHB= (R1+TP*SVP1) /BK 
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ZHB= ( Z 1+TP*CVP1 ) /BK 
TH=ATAN2 ( RHB , ZHB+ 1 . E- 5 ) 

CC=COS (Vl-TH) 

SS=SIN (Vl-TH) 

RR=SQRT (RHB**2+ZHB**2 ) 

S 1=S 1+AT ( L) * ( . 5+TP/D1 ) *CC*CIRCRHO ( CNPHI , XP , AP , CNRHO , XR , AR , 

*  ARAD, SCAN, PHI , RHB, ZHB) 

S2=S2+AT (L) * ( . 5+TP/D1) *SS*CIRCTHETA (CNPHI , XP, AP, CNRHO , XR, AR , 

*  ARAD, SCAN, PHI, RHB,ZHB) 

13  CONTINUE 
Sl=Sl*Dl/2. 

S2=S2*Dl/2. 

C  t-CURRENT  INTEGRATION  FOR  THE  NEGATIVE  SLOPE 
DO  14  L= 1 , NT 
TP=D2  *XT ( L) / 2 . 

RHB= (R2+TP*SVP2 ) /BK 
ZHB= ( Z2+TP*CVP2 ) /BK 
TH=ATAN2 (RHB, ZHB+1 . E-5) 

SS=SIN (V2-TH) 

CC=COS (V2-TH) 

RR=SQRT(RHB**2+ZHB**2) 

S  3 =S  3 + AT ( L) * ( . 5 -TP/ D2 ) *  CC*CIRCRHO ( CNPHI , XP , AP , CNRHO , XR , AR , 

*  ARAD, SCAN, PHI, RHB, ZHB) 

S4=S4+AT (L) * ( . 5-TP/D2 ) *SS*CIRCTHETA (CNPHI , XP, AP, CNRHO, XR, AR, 

*  ARAD, SCAN, PHI, RHB, ZHB) 

14  CONTINUE 
S3=S3*D2/2 . 

S4=S4  *D2/2 . 

ENDIF 

C  phi-CURRENT  INTEGRATION 
DO  15  L=1 , NT 
TP=D1*XT(L) /2. 

RHB=(R1+TP*SVP1)/BK 
ZHB=(Z1+TP*CVP1)/BK 
TH=ATAN2 (RHB , ZHB+1 . E-5 ) 

RR=SQRT (RHB**2+ZHB**2 ) 

S5=S5+AT ( L) *CIRCPHI (CNPHI , XP, AP , CNRHO , XR , AR , 

*  ARAD, SCAN , PHI , RHB, ZHB) /Rl* (RHB*BK) 

15  CONTINUE 
S5=S5*Dl/2 . 

SP=SP+A(J) *PSI*S5 
IF (IP. LE.MT)  THEN 
ST1=ST1+A(J) *PSI*S1 
ST2=ST2+A ( J) *PSI*S2 
ST3=ST3+A(J) *PSI*S3 
ST4=ST4+A ( J) *PSI*S4 
ENDIF 
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20  CONTINUE 

C  COMPONENTS  ARE  STORED  IN  A  COLUMN  VECTOR:  VT(1,MT), 
C  VP(MT+1,N) 

IF(IP.LE.MT)  THEN 

R(IP) = (ST1+ST2+ST3+ST4 ) *P1 

ENDIF 

R(IP+MT)=SP*P1 
30  CONTINUE 

RETURN 
END 


C  *  *  *  *  *  *  *FUNCTION  CIRCTHETA. 

C  FUNCTION  :  CIRCTHETA 

C  DATE  :  1  OCTOBER  1991 

C  REVISED  :  17  February  1992 

C  PROGRAMMER  :  R.M.  FRANCIS 

C  COMMENTS  :  THIS  FUNCTION  WILL  RIGOROUSLY  CALCULATE  THE 
C  ELECTRIC  FIELD  FOR  A  CIRCULAR  APERTURE.  THE  FIELD  IS 
C  CALCULATED  AT  THE  BOUNDARY  DEFINED  BY  AN  OGIVE.  THE  APERTURE 
C  IS  LOCATED  AT  Z  =  0,  AND  IS  BORE  SIGHTED  ON  THE  Z  AXIS.  THE 
C  SUBROUTINE  OGIVE  IS  THE  SOURCE  OF  THE  GEOMETRIC  DATA 
C  REQUIRED  BY  CIRC  TO  PERFORM  COMPUTATIONS. 

C  ALL  PHYSICAL  DIMENSIONS  ARE  NORMALIZED  TO  WAVELENGTH. 


COMPLEX  FUNCTION  CIRCTHETA  ( CNPHI ,  XP ,  AP ,  CNRHO ,  XR ,  AR ,  ARAD ,  SCAN , 

*  PHI , RHB, ZHB) 

INTEGER  CNPHI, CNRHO, M,N 

REAL  XP(100) , AP(100) ,XR(50) ,AR(50) , 

*  ARAD, SCAN, PHI, RHB, ZHB 
C  LOCAL  VARIABLES 

REAL  Z , RZ , PHIPRIME , PI , S , C, 

*  R,RP,X2, Yl, Y2,K 
COMPLEX  J , JK, G1 , G2 , XI , CON, CC, 

*  sumx , sumy , sumz , CIRCTHETA 
RR1= (ARAD+0 . ) /2 . 

RR2= (ARAD-0. )/2. 

PI=3. 141592654 
K=6. 283185307 
J  =  (0.0, 1.0) 

JK=  (0.0,6.283185307) 

CON= (0.0, -15.0) 

CC=  C0N*RR1 

C  OUTER  LOOP:  INTEGRATE  OVER  PHI .  .  . -PKPHKPI . 

C  INNER  LOOP:  INTEGRATE  OVER  RHO...  0  <RHO<ANT 
RZ=RHB 
Z=ZHB 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 
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sumx= ( 0 . 0 , 0 . 0 ) 
sumy= ( 0 . 0 , 0 . 0 ) 
sumz=(0.0,0.0) 

DO  50  M=1 , CNPHI 

PHIPRIME=PI*XP(M) 

DO  60  N=1 , CNRHO 
RP=RR1*XR(N)+RR2 

R=SQRT ( ( (RZ-RP)**2+Z**2)+(4*RZ*RP*SIN( ( PHI-PHIPRIME) /2 ) **2) ) 
Gl= ( ( (K*R)**2)-1-(JK*R) )/R**3 
G2=(3+(3*JK*R)-(K*R) **2)/R**5 
X1=JK* ( (RP*COS (PHIPRIME) *SIN (SCAN) ) -R) 

Yl= (RZ*COS (PHI) )-(RP*COS (PHIPRIME) ) 

X2=Y1**2 

Y2= (RZ*SIN (PHI) ) -(RP*SIN (PHIPRIME) ) 

sumx=sumx+ (CC*AP(M) *AR(N) * (G1+(X2*G2) ) *CEXP (XI ) *RP) 
sumy=sumy+ (CC*AP(M) *AR(N) *Y1*Y2*G2*CEXP(X1) *RP) 
sumz=sumz+(CC*AP(M) *AR (N) *Y1*Z*G2*CEXP(X1) *RP) 

60  CONTINUE 
50  CONTINUE 

CIRCTHETA= ( C* COS ( PHI ) *  sumx )  +  ( C*S IN ( PHI ) *  sumy ) - ( S  *  sum z ) 

RETURN 

END 


C***********FUNCTION  CIRCRHO. 

C  FUNCTION  :  CIRCRHO 

C  DATE  :  1  OCTOBER  1991 

C  REVISED  :  17  February  1992 

C  PROGRAMMER  :  R.M.  FRANCIS 


COMPLEX  FUNCTION  CIRCRHO (CNPHI , XP, AP, CNRHO, XR, AR, ARAD, 

*  SCAN, PHI, RHB,ZHB) 

INTEGER  CNPHI, CNRHO, M,N 

REAL  XP(100) , AP ( 100) ,XR(50) ,AR(50) , 

*  ARAD, SCAN, PHI, RHB,ZHB 
C  LOCAL  VARIABLES 

REAL  Z,RZ, PHIPRIME, PI, S,C, 

*  R,RP,X2,Y1,Y2,K 
COMPLEX  J, JK,G1,G2,X1,C0N,CC, 

*  sumx, sumy , sumz , CIRCRHO 
RR1= (ARAD-0. )/2. 

RR2=(ARAD+0. )/2. 

PI=3. 141592654 
K=6. 283185307 
J  =  (0.0, 1.0) 

JK=  (0.0,6.283185307) 

CON=(0. 0, -15. 0) 

CC=C0N*RR1 

C  OUTER  LOOP:  INTEGRATE  OVER  PHI .  .  .  -PKPHKPI . 


C  INNER  LOOP:  INTEGRATE  OVER  RHO. . .  0  <RHO<ANT 
RZ=RHB 
Z=ZHB 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 
sumx=(0.0,0.0) 
sumy=(0.0,0.0) 
sumz=(0. 0,0.0) 

DO  50  M=1 , CNPHI 

PHIPRIME=PI *XP (M) 

DO  60  N=1 , CNRHO 
RP=RR1*XR(N)+RR2 

R=SQRT ( ( (RZ-RP) **2+Z**2) + (4*RZ*RP*SIN( (PHI-PHIPRIME) /2 ) **2 ) ) 
Gl= ( ( (K*R) **2)-l-(JK*R) )/R**3 
G2= ( 3+ ( 3  *JK*R) - (K*R) *  *2 ) /R**5 
X1=JK* ( (RP*COS (PHIPRIME) *SIN (SCAN) ) -R) 

Yl= (RZ*COS (PHI) ) - (RP*COS (PHIPRIME) ) 

X2=Y1**2 

Y2= (RZ*SIN (PHI) )-(RP*SIN (PHIPRIME) ) 

svnnx=sumx+ ( CC*AP (M) *AR(N) * (G1+ (X2*G2) ) *CEXP(X1) *RP) 
sumy=sumy+(CC*AP(M) *AR(N) *Y1*Y2*G2*CEXP(X1) *RP) 
sumz=sumz+(CC*AP(M) *AR(N) *Y1*Z*G2*CEXP(X1) *RP) 

60  CONTINUE 
50  CONTINUE 

CIRCRHO= ( S*COS ( PHI ) *sumx) + (S*SIN (PHI ) *sumy ) + ( C*sumz ) 

RETURN 

END 


C*******FUNCTION  CIRCPHI. 

C  FUNCTION  :  CIRCPHI 

C  DATE  :  1  OCTOBER  1991 

C  REVISED  :  17  Feb  1992 

C  PROGRAMMER  :  R.M.  FRANCIS 

COMPLEX  FUNCTION  CIRCPHI (CNPHI, XP,AP, CNRHO, XR,AR, ARAD, 

*  SCAN, PHI, RHB,ZHB) 

INTEGER  CNPHI, CNRHO, M,N 

REAL  XP(100) , AP (100) ,XR(50) ,AR(50) , 

*  ARAD, SCAN, PHI, RHB,ZHB 
C  LOCAL  VARIABLES 

REAL  Z,RZ, PHIPRIME, PI, S,C, 

*  R, RP, X2 , Y1 , Y2 , K, RR1 , RR2 
COMPLEX  J, JK,G1,G2,X1,C0N,CC, 

*  sumx, sumy, sumz 
PI=3. 141592654 

K=6 . 283185307 
RR1= (ARAD-0. )/2. 

RR2= ( ARAD+0 . ) /2 . 

J  =  (0.0, 1.0) 

JK=  (0.0,6.283185307) 
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CON=(0.0,-15.0) 

CC=C0N*RR1 

C  OUTER  LOOP:  INTEGRATE  OVER  PHI .  .  .  -PKPHKPI . 

C  INNER  LOOP:  INTEGRATE  OVER  RHO. . .  0  <RHO<ANT 
RZ=RHB 
Z=ZHB 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 
surox=(0.0,0.0) 
sumy= ( 0 . 0 , 0 . 0 ) 
sumz=(0. 0, 0. 0) 

DO  50  M=1 , CNPHI 

PHIPRIME=PI*XP(M) 

DO  60  N=1 , CNRHO 
RP=RR1*XR (N) +RR2 

R=SQRT ( ( (RZ-RP) **2+Z**2 ) + (4*RZ*RP*SIN ( (PHI -PHI PRIME) /2 ) **2 ) ) 
Gl=( ( (K*R) **2)-l-(JK*R) )/R**3 
G2= (3+ (3*JK*R) - (K*R) **2)/R**5 
X1=JK* ( (RP*COS (PHIPRIME) *SIN(SCAN) ) -R) 

Yl=(RZ*COS (PHI) ) - (RP*COS (PHIPRIME) ) 

X2=Y1**2 

Y2= (RZ*SIN (PHI) ) - (RP*SIN (PHIPRIME) ) 

sumx=sumx+ ( CC*AP (M) *AR(N) * (G1+ (X2*G2) ) *CEXP(X1) *RP) 
sumy=sumy+(CC*AP(M) *AR(N) *Y1*Y2*G2*CEXP (XI) *RP) 
sumz=sumz+ ( CC*AP (M) *AR(N) *Y1*Z*G2*CEXP(X1) *RP) 

60  CONTINUE 
50  CONTINUE 

CIRCPHI= (-SIN (PHI) *sumx) + (COS (PHI) *sumy) 

RETURN 

END 


C*********SUBROUTINE  TESTSPHERE 


C 

C 

C 

C 

C 


SUBROUTINE 

DATE 

PROGRAMMER 

COMMENTS 


TESTSPHERE 
21  FEBRUARY  1992 
R.  M.  FRANCIS 

GENERATES  A  SPHERE,  WITH  PLOTTED  POINTS 
AT  EQUAL  INTERVALS  IN  THETA. 


SUBROUTINE  TESTSPHERE (NP, ZH , RH , BASE , RS , ZP) 
INTEGER  NP, I 

REAL  ZH(400) ,RH(400) , BASE , RS , ZP, SPRAD 
PI=3. 1435926 
BK=2 . *PI 
ZP=0.0 

WRITE (6,*) 'ENTER  AN  EVEN  NUMBER  OF  POINTS  (NP) : ' 
READ ( 5 , *) NP 

WRITE (6, *) 'ENTER  SPHERE  RADIUS' 

READ (5,*) SPRAD 
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BASE=SPRAD 
RS=SPRAD 
DO  1241  1=1, NP 

ANGLE=PI* FLOAT ( 1-1) / ( 2 . * FLOAT (NP-1) ) 

ZH (I) =BK*SPRAD*COS (ANGLE) 

RH ( I) =BK*SPRAD*SIN ( ANGLE) 

1241  CONTINUE 
RETURN 
END 

B.  SAMPLE  DATA  FILE  OPTLDBOR 

***  BOR  RADIATION  PATTERN  FOR  A  CIRCULAR  DISC  RADIATOR  USING 
GENEX  *** 


BOR  DIAMETER  (WAVELENGTHS) =***** 
NUMBER  OF  GENERATING  POINTS  (NP)=  70 

SURFACE  RADIUS50 . 00  2PRIME  0.00 


INDEX 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 
27 


NT  NPHI 
4  60 


Z(I) 

RHO(I) 

50.000 

0.000 

49.987 

1.138 

49.948 

2.276 

49.883 

3.412 

49.793 

4.547 

49.676 

5.679 

49.534 

6.808 

49.366 

7.934 

49.173 

9.056 

48.954 

10.173 

48.710 

11.285 

48.440 

12.390 

48.146 

13.490 

47.826 

14.582 

47.482 

15.667 

47.113 

16.744 

46.720 

17.812 

46.302 

18.871 

45.861 

19.920 

45.395 

20.959 

44.906 

21.987 

44.394 

23.003 

43.859 

24.008 

43.301 

25.000 

42.721 

25.979 

42.119 

26.945 

41.494 

27.897 

SURF 

IMPED 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 
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28 

40.848 

28.834 

10000.00 

0.00 

29 

40.182 

29.756 

10000.00 

0.00 

30 

39.494 

30.663 

10000.00 

0.00 

31 

38.786 

31.554 

10000.00 

0.00 

32 

38.057 

32.429 

10000.00 

0.00 

33 

37.309 

33.287 

10000.00 

0.00 

34 

36.542 

34.128 

10000.00 

0.00 

35 

35.755 

34.951 

10000.00 

0.00 

36 

34.951 

35.755 

10000.00 

0.00 

37 

34.128 

36.542 

10000.00 

0.00 

38 

33.287 

37.309 

10000.00 

0.00 

39 

32.429 

38.057 

10000.00 

0.00 

40 

31.554 

38.786 

10000.00 

0.00 

41 

30.663 

39.494 

10000.00 

0.00 

42 

29.756 

40.182 

10000.00 

0.00 

43 

28.834 

40.848 

10000.00 

0.00 

44 

27.897 

41.494 

10000.00 

0.00 

45 

26.945 

42.119 

10000.00 

0.00 

46 

25.979 

42.721 

10000.00 

0.00 

47 

25.000 

43.301 

10000.00 

0.00 

48 

24.008 

43.859 

10000.00 

0.00 

49 

23.003 

44.394 

10000.00 

0.00 

50 

21.987 

44.906 

10000.00 

0.00 

51 

20.959 

45.395 

10000.00 

0.00 

52 

19.920 

45.861 

10000.00 

0.00 

53 

18.871 

46.302 

10000.00 

0.00 

54 

17.812 

46.720 

10000.00 

0.00 

55 

16.744 

47.113 

10000.00 

0.00 

56 

15.667 

47.482 

10000.00 

0.00 

57 

14.582 

47.826 

10000.00 

0.00 

58 

13.490 

48.146 

10000.00 

0.00 

59 

12.390 

48.440 

10000.00 

0.00 

60 

11.285 

48.710 

10000.00 

0.00 

61 

10.173 

48.954 

10000.00 

0.00 

62 

9.056 

49.173 

10000.00 

0.00 

63 

7.934 

49.366 

10000.00 

0.00 

64 

6.808 

49.534 

10000.00 

0.00 

65 

5.679 

49.676 

10000.00 

0.00 

66 

4.547 

49.793 

10000.00 

0.00 

67 

3.412 

49.883 

10000.00 

0.00 

68 

2.276 

49.948 

10000.00 

0.00 

69 

1.138 

49.987 

10000.00 

0.00 

70 

0.000 

50.000 

10000.00 

0.00 

PHI  OF  RECEIVER  (DEG) =  0.00 

MAXIMUM  FIELD  VALUE  (V/M)  =  0.16659E+03 


98 


ANGLE 

CO- 

■POL 

X-POL 

(DEG) 

(VOLTS) 

(DEG) 

(DB-REL) 

(VOLTS) 

(DEG) 

(DB-REL) 

0.00 

6.9271 

-104.51 

-27.62 

0.0000 

-92.95 

-100.00 

1.30 

3.6864 

-106.76 

-33.10 

1.1253 

131.93 

-43.41 

2.61 

2.6516 

-50.40 

-35.96 

0.7957 

163.87 

-46.42 

3.91 

4.8350 

-14.07 

-30.74 

0.5577 

-172.05 

-49.50 

5.22 

6.0401 

-5.71 

-28.81 

0.2076 

-139.66 

-58.09 

6.52 

4.7465 

-2.93 

-30.91 

0.2741 

5.44 

-55.68 

7.83 

1.3824 

-4.09 

-41.62 

0.5842 

30.18 

-49.10 

9.13 

2.8075 

-174.09 

-35.47 

0.7556 

57.97 

-46.87 

10.43 

5.8878 

-174.09 

-29.03 

0.9056 

82.84 

-45.29 

11.74 

6.3663 

-172.38 

-28.36 

0.9042 

101.87 

-45.31 

13.04 

3.9668 

-170.31 

-32.46 

0.6473 

118.37 

-48.21 

14.35 

0.9093 

7.83 

-45.26 

0.2436 

-130.52 

-56.70 

15.65 

6.0107 

6.17 

-28.85 

1.0247 

-96.02 

-44.22 

16.96 

8.2390 

4.51 

-26.12 

1.2154 

-95.19 

-42.74 

18.26 

6.7229 

7.30 

-27.88 

0.9426 

-80.30 

-44.95 

19.57 

2.3749 

15.98 

-36.92 

0.4040 

-54.39 

-52 .31 

20.87 

3.7384 

-170.33 

-32.98 

0.5796 

72.59 

-49.17 

22.17 

9.6622 

-170.82 

-24.73 

1.3618 

88.39 

-41.75 

23.48 

10.7713 

-175.50 

-23.79 

1.5721 

81.57 

-40.50 

24.78 

6.6302 

-160.18 

-28.00 

0.6277 

115.25 

-48.48 

26.09 

3.5922 

-149.07 

-33.33 

0.8164 

137.53 

-46.20 

27.39 

7.2792 

27.16 

-27.19 

0.9508 

-75.15 

-44.87 

28.70 

13.9956 

3.03 

-21.51 

2.1212 

-107.94 

-37.90 

30.00 

12.6783 

26.06 

-22.37 

1.1997 

-71.27 

-42.85 

31.30 

14.4138 

12.62 

-21.26 

2.0480 

-86.02 

-38.21 

32.61 

1.4466 

93.02 

-41.23 

0.5246 

44.99 

-50.04 

33.91 

5.0101 

-168.15 

-30.44 

0.7135 

42.79 

-47.37 

35.22 

18.2410 

-161.85 

-19.21 

2.3460 

82.73 

-37.03 

36.52 

22.8067 

-158.28 

-17.27 

2.5909 

93.03 

-36.16 

37.83 

17.8045 

-160.35 

-19.42 

2.0096 

80.80 

-38.37 

39.13 

18.0824 

-146.48 

-19.29 

2.4138 

124.03 

-36.78 

40.43 

9.1689 

3.14 

-25.19 

1.7552 

-103.63 

-39.55 

41.74 

13.8129 

38.91 

-21.63 

0.9059 

-103.47 

-45.29 

43.04 

30.3979 

9.45 

-14.78 

4.3773 

-120.67 

-31.61 

44 .35 

47.4641 

29.97 

-10.91 

5.2557 

-84.81 

-30.02 

45.65 

36.1417 

33.25 

-13.27 

3.3511 

-87.55 

-33.93 

46.96 

31.6252 

32.13 

-14.43 

3.2704 

-81.37 

-34 . 14 

48.26 

14.5891 

62.48 

-21.15 

2.0676 

-25.23 

-38.12 

49.57 

18.9270 

-170.31 

-18.89 

2.8516 

41.68 

-35.33 

50.87 

45.8795 

-157.68 

-11.20 

5.2496 

54.10 

-30.03 

52.17 

78.0162 

-154.23 

-6.59 

8.1824 

59.69 

-26.18 

53.48 

106.099 

-152.60 

-3.92 

10.9266 

61.19 

-23.66 

54.78 

136.447 

-152.71 

-1.73 

13.7615 

65.23 

-21.66 

56.09 

160.822 

-147.27 

-0.31 

14.3460 

71.03 

-21.30 

57.39 

161.358 

-145.24 

-0.28 

13.9190 

68.77 

-21.56 

58.70 

166.587 

-146.15 

0.00 

14.2015 

67.31 

-21.39 

60.00 

162.565 

-144.37 

-0.21 

13.4435 

67.27 

-21.86 

61.30 

153.157 

-144.83 

-0.73 

12.5463 

66.68 

-22.46 

62.61 

143.692 

-144.21 

-1.28 

11.3471 

67.74 

-23 . 34 
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65.22 

108.369 

-140.84 

-3.73 

7.6142 

65.73 

-26.80 

66.52 

85.0397 

-139.06 

-5.84 

5.6224 

60.07 

-29.43 

67.83 

67.7272 

-140.66 

-7.82 

4.2600 

58.69 

-31.84 

69.13 

59.7937 

-143.81 

-8.90 

3.8263 

67.63 

-32.78 

70.43 

47.3499 

-145.20 

-10.93 

3.1370 

62.89 

-34 . 50 

71.74 

30.9244 

-145.07 

-14.63 

2.1385 

43.54 

-37.83 

73.04 

19.4698 

-143.32 

-18.65 

1.2247 

37.25 

-42.67 

74.35 

12.7001 

-149.46 

-22.36 

0.8089 

32.28 

-46.28 

75.65 

5.2649 

-153.77 

-30.01 

0.8572 

-15.23 

-45.77 

76.96 

8.7952 

-119.58 

-25.55 

0.7979 

-159.71 

-46.39 

78.26 

12.4865 

173.28 

-22.50 

1.4228 

77.84 

-41.37 

79.57 

11.0782 

91.40 

-23 . 54 

1.4678 

-11.88 

-41.10 

80.87 

3.7571 

-1.39 

-32  .  < 

1.0531 

-129.37 

-43.98 

82.17 

6.6444 

103.97 

-27.9c 

0604 

30.52 

-43.92 

83.48 

6.3822 

-4.11 

-28.33 

68 

-128.96 

-43.87 

84.78 

3.0249 

-158.22 

-34.82 

'4 

104.94 

-54.94 

86.09 

1.8402 

75.64 

-39.14 

0.2566 

-5.90 

-56.25 

87.39 

3.9405 

-150.16 

-32.52 

0.8277 

166.07 

-46.08 

88.70 

7.3638 

110.06 

-27.09 

1.1926 

21.33 

-42.90 

90.00 

7.6190 

-5.04 

-26.79 

1.2270 

-108.00 

-42.66 

C.  SOURCE  CODE  FOR  GENPLT.M 

% PROGRAM  :  genplt.m 

%DATE  :  15  February  1992 

% PROGRAMMER:  R.M.  FRANCIS 

%Comments  :  This  prc-ram  is  ?  rlottirr  routine  for  the  data 
%  points  generated  :  -  •  -  -rv  rrar  T.D3GRMM.F  and  the 

%  program  TEST.M.  tne  runs  of  test.m  and 

%  LDBORMM . F . 


load  bang.m 
load  bxpole.m 
load  bcpole.m 
load  ctest 
for  i=l:np 

if  bcpole ( i) <=0 . 0 
bcpole ( i) =eps ; 

end 

if  bxpole ( i) <=0 . 0 
bxpole (i)=eps; 

end 

el (i) =10 *log (bcpole (i) )  ; 
ep(i)=10*log (bxpole (i) ) ; 
zl (i)=10*log(z (i) )  ; 
if  zl(i)<-50 
zl (i) =-50; 

end 
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if  el(i)<-50 
el (i)=-50; 

end 

if  ep(i)<-50 
ep(i)=-50; 

end 

end 

subplot(221) ,plot(deg, z , bang, bcoole, ' . ') 
title ( 'TEST  OF  LDBORMM .  F ' )  , xlabel ( ' degrees ' ) 
ylabel ( 'magnitude' ) 

subplot  (222)  ,  plot  (deg,  zl,b_ng,el,  ' . ') 
title ('LOG  PLOT  (db) ') , xlabel ( 'degrees' ) 
ylabel ( 'magnitude' ) 

text (. 1 ,. 4 , 'analytical  (solid),  LDBORMM. F( .)', 'sc' ) 
text ( . 1 , . 3 , ' PARAMETERS : ' , ' sc ' ) 
text (. 1 ,. 25 , 'number  of  points  =  90', 'sc') 
text ( .5, .25, 'scan  angle  =  60  deg.', 'sc') 
text ( . 5 , . 15 , ' CNRHO  =  20', 'sc') 
text (.5, .2, 'CNPHI  =  60', 'sc') 
text ( . 1 , . 2 , 'NPHI  =  100', 'sc') 
text(.l, .15, 'NT  =  10', 'sc') 
text (. 1, . 1, 'distance  =  50', 'sc') 
text (. 5 ,. 1 , 'aperture  radius  =  5', 'sc') 
text (. 1 ,. 05 , 'complex  impedance  =  (10000.0,  0.0)', 'sc') 
pause 
%print 
clg 

load  etscat.m 
load  epscat.m 
load  etf.m 
load  epf.m 

subplot (221) , plot (bang, etscat) 
title ( 'ETSCAT') 

subplot(222) , plot (bang, epscat) 
title ( 'EPSCAT' ) 
subplot(223) , plot (bang, etf) 
title ( 'ETF' ) 

subplot(224) , plot (bang, epf) 
title ( 'EPF' ) 

%print 
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APPENDIX  B 

TABLE  B.l.  VARIABLES  IN  SUBROUTINE  GENEXX 


VARIABLE 

TYPE 

DESCRIPTION 

RANGE 

NM 

integer 

Number  of  modes  in  0 

1<  NM  <  oo 

NP 

integer 

Even  number  of  points  on 
surface  of  the  BOR 

2<  NP  <400 

NT 

integer 

Even  number  of  points  of 
integration  in  T  on  BOR 

2<  NT  <4 

NPHI 

integer 

Even  number  of  points  of 
integration  in  0  on  BOR 

2<NPHI<100 

CNRHO 

integer 

Even  number  of  points  of 
integration  ir.  p' (source) 

2<CNRHO<50 

CNPHI 

integer 

Even  number  of  points  of 
integration  in  0*  (source) 

2<CNPHI<100 

XT,X,XR, 

XP 

real 

Abscissa  of  the  function  (xn) 
in  T,0,p'  and  0' 

-i<xn<i 

AT , A , AR , 
AP 

real 

Weights  of  the  series 
terms (for  xQ)  approximating 
the  integral 

0<t»)n<l 

SCAN 

real 

angular  displacement  (in 
degrees)  of  source  boresight 
and  Z  axis 

-9O<0s<9O 

PHI 

real 

orientation  in  <p  of  receiver 

O<0<18O 

ARAD 

real 

dimensionless  radius  of 
source  antenna  (a/A) 

0<p  '  <  oo 

RH(i) 

real 

the  product  of  the  wave 
number  (k)  and  the  radius  of 
the  BOR,  units  in  radians 

0<  r  <  oo 

3<  i  <  399 

ZH  ( i ) 

real 

the  product  of  the  wave 
number  (k)  and  the  distance  z 
on  the  BOR  axis,  units  in 
radians 

0<  z  <  oo 

3<  i  <  399 

R(i) 

complex 

the  excitation  vector 
i  =2NP— 3 

max 

-  oo<  R(i)  <oo 

4 <  i  <  1600 
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APPENDIX  C 

A.  PROGRAM  TO  GENERATE  DATA  FILES  FOR  GAUSS  QUADRATURE 
INTEGRATION 

C  PROGRAM  :  GAUS.F 

C  DATE  :  12  SEPTEMBER  1991 

C  PROGRAMMER  :  R.  FRANCIS,  D.  JENN 

C  PROGRAM  TO  COMPUTE  WEIGHTS  AND  ABSCISSAS  FOR  GAUSSIAN  INT. 
CHARACTER*?  FNAME 
DIMENSION  X ( 500) , A (500 ) 

WRITE (6 , *)  'ENTER  FILE  NAME  (gaus - )' 

READ (5,*) FNAME 

WRITE (6 , *)  'ENTER  AN  EVEN  NUMBER  (max  pts  is  500)' 

READ ( 5 , * )  N 

OPEN ( 2 , FILE=FNAME) 

WRITE (2 , * )  N 

CALL  GAUSLEG (-1. ,1. ,X,A,N) 

DO  10  1=1, N 

10  WRITE ( 2 , * )  X(I) ,A(I) 

STOP 

END 

SUBROUTINE  GAUSLEG (XI , X2 , X, W, N) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

REAL* 4  X1,X2,X(N) ,W(N) 

PARAMETER  (EPS=3.D-14) 

M=(N+l)/2. 

XM=. 5D0* (X2+X1) 

XL= . 5D0* (X2-X1 ) 

DO  12  1=1, M 

Z=COS (3.1415926535D0*(I~. 25D0) / (N+ . 5 DO) ) 

1  CONTINUE 

Pl=l . DO 
P2=0 . DO 
DO  11  J=1 , N 
P3=P2 
P2=P1 

Pl=( (2.D0*J-1.D0) *Z*P2- ( J-l . DO) *P3)/J 

11  CONTINUE 

PP=N* (Z*P1-P2)/ (Z*Z-1. DO) 

Z1=Z 

Z=Z1-P1/PP 

IF(ABS(Z-Z1) .GT.EPS)  GO  TO  1 
X ( I ) =XM-XL*Z 
X (N+l-I ) =XM+XL*Z 
W(I)=2.D0*XL/( (l.D0-Z*Z) *PP*PP) 

W (N+l-I ) =W ( I ) 

12  CONTINUE 
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RETURN 

END 


B.  SAMPLE  DATA  FILE  GENERATED  BY  GAPS . F 

The  following  is  an  example  of  the  external  sequential 
data  file  generated  by  the  program  GAUS.F.  The  first  row  is  an 
even  integer  (N)  .  The  data  in  rows  two  through  N  are  the 
abscissa  and  weights  for  the  series  approximation  of  the 
integral  calculation. 

Sample  data  file;  qaus20 _ 


20 

-0.993129 

1 . 76140E-02 

-0.963972 

4 . 06014E-02 

-0.912234 

6 . 26720E-02 

-0.839117 

8 . 32767E-02 

-0.746332 

1. 01930E-01 

-0.636054 

0.118195 

-0.510867 

0.131689 

-0.373706 

0.142096 

-0.227786 

0.149173 

-7 . 65265E-02 

0.152753 

7 . 65265E-02 

0.152753 

0.227786 

0.149173 

0.373706 

0.142096 

0.510867 

0.131689 

0.636054 

0.118195 

0.746332 

1 . 01930E-01 

0.839117 

8 . 32767E-02 

0.912234 

6 . 26720E-02 

0.963972 

4 . 06014E-02 

0.993129 

1 . 76140E-02 
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APPENDIX  D 


A.  ARGUMENTS:  CIRCTHETA.  CIRCRHO  AND  CIRCPHI 
TABLE  D.l.  VARIABLE  LIST  FOR  FUNCTION  SUBPROGRAMS . 


VARIABLE 

TYPE 

DESCRIPTION 

RANGE 

CNPHI 

integer 

see  Table  B.l. 

same 

CNRHO 

integer 

see  Table  B.l 

same 

XP,  XR 

real 

see  Table  B.l 

same 

AP,  AR 

real 

see  Table  B. 1 

same 

ARAD 

real 

see  Table  B.l 

same 

SCAN 

real 

see  Table  B.l 

same 

PHI 

real 

see  Table  B.l 

same 

RHB 

real 

same  as  RH(i)  of 

Table  B.l 

same 

ZHB 

real 

same  as  ZH(i)  of 

Table  B.l 

same 

circtheta 

circrho 

circphi 

complex 

Eg,  Ep  and  Ep 
calculated  at  the 
point  on  the  BOR 
surface  defined  by 

RHB, ZHB  and  PHI 

all 

complex 

numbers 

B.  TEST  PROGRAM:  CIRCSUB.F 

The  following  program  was  used  to  test  for  convergence  the 
numerical  methods  used  to  calculate  the  electric  field 
components  for  a  uniformly  illuminated  circular  aperture. 
Before  executing  the  main  program  LDBORMM. F,  it  is  desirable 
to  determine  the  minimum  number  of  terms  required  (CNPHI  and 
CNRHO)  for  the  series  approximation  of  the  source  field  to 
converge  to  the  analytical  solution  of  Equation  (3.26).  The 
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following  programs  are  provided  to  test  for  the  number  of 
terms  required  for  convergence. 


The  programs  CIRCSUB.F,  TEST.M  and  SCANPLT.M  are  executed 
in  the  same  directory.  To  minimize  the  number  of  files  in  the 
main  directory,  maintain  these  programs  in  a  subdirectory 
named  "circtest" .  CIRCSUB.F  creates  the  external  sequential 
data  files  ang.m,  etheta.m,  erho.m  and  ephi.m.  These  data 
files  are  vectors  of  length  NP,  containing  the  angle  (0)  and 
the  magnitude  of  E0,  Ep  and  E^.  The  matlab  program  TEST.M 
calculates  the  analytic  solution  for  Efl  in  accordance  Equation 
(3.26)  and  writes  the  data  to  a  file  CTEST . MAT .  The  matlab 
program  SCANPLT.M  loads  the  data  files  and  plots  the  results. 

1.  Source  Code  for  CIRCSUB.F 


C  TESTPROGRAM  :  CIRCSUB.F 

C  DATE  :  1  OCTOBER  1991 

C  REVISED  :  12  JANUARY  1992 

C  PROGRAMMER  :  R.M.  FRANCIS 

C  COMMENTS  :  THIS  SUBROUTINE  CALCULATES  THE  ELECTRIC  FIELD 
C  FOR  A  CIRCULAR  APERTURE.  THE  FIELD  IS  CALCULATED  AT  THE 
C  BOUNDARY  DEFINED  BY  AN  OGIVE.  THE  APERTURE  IS  LOCATED  AT  Z 
C  =  0,  AND  IS  BORE  SIGHTED  ON  THE  Z  AXIS.  THE  SUBROUTINE 
C  OGIVE  IS  THE  SOURCE  OF  THE  GEOMETRIC  DATA  REQUIRED  BY  CIRC 
C  TO  PERFORM  COMPUTATIONS. 

C  ALL  PHYSICAL  DIMENSIONS  ARE  NORMALIZED  TO  WAVELENGTH. 


C*********  PROGRAM  CIRC 

C  DECLARE  VARIABLES  INTEGER  NP , NPHI , NRHO , I , M , N 
C  THE  CHARACTER  STRING  PHIPTS/RHOPTS  WILL  ACCESS  THE 
C  APPROPIATE  FILE  WITH  THE  X'S  AND  COEEFFICIENTS ,  TO  PERFORM 
C  INTEGRATION  BY  GAUSSIAN  QUADRATURE.  EG.  GAUS10  YEILDS ,  10 
C  POINTS ;GAUS20,  20PTS ,  ETC. 

CHARACTER* 15  PHIPTS , RHOPTS , GP, GR 

REAL  XP(400) ,AP(400) ,XR(100) ,AR(100) , S , C, P, PI , SPRAD, 

*  SCAN, Z,RZ, PHI, PHIPRIME, DEG, SA, K, ET, ER, EP, 
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*  R,  RP, X2 , Y1 , Y2 , ARAD, RH ( 4  00) ,ZH(400) ,angle(400) 
COMPLEX  J, JK,G1,G2,X1,C0N,CC, 

*  sumx , sumy , sumz , 

*  ETHETA(400) ,ERHO(400) ,EPHI(400) 

PI=3. 141592654 

J  =  (0.0, 1.0) 

JK=  (0.0,6.283185307) 

K=6. 283185307 
CON= (0.0, -15 . 0) 

WRITE (6,*) 'ENTER  THE  RADIUS  OF  THE  ANTENNA, (in  wavelengths) ' 
READ ( 5 , * )  ARAD 
ANT=ARAD/2 . 

WRITE ( 6 ,*)• ENTER  PHI , (OBSERVATION)  IN  DEGREES' 

READ ( 5 , * ) P 
PHI=P*PI/180 . 

WRITE (6,*) 'ENTER  SCAN  ANGLE  IN  DEGREES' 

READ ( 5 , * ) SA 
SCAN=-SA*PI/180. 

WRITE(6,*) 'ENTER  PHIPTS  STRING:  "gaus ###  " , (max  400)' 
READ ( 5 , * ) GP 
PHIPTS= ' gaus/ ' //GP 

WRITE(6,*) 'ENTER  RHOPTS  STRING:  "gaus###  ", (max  100)' 
READ ( 5 , * ) GR 
RHOPTS= ' gaus/ • //GR 
CALL  OGIVE (NP,ZH,RH) 
make  loop  to  enter  test  sphere. 

write (6 ,*)' enter  the  sphere  radius' 
read ( 5 , * ) SPRAD 

write ( 6 ,*)' enter  the  number  of  points' 
read (5 , *) NP 
do  888  1=1, NP 

angle (I) =PI* (1-1) / (2* (NP-1) ) 

RH (I) =SPRAD*SIN (angle (I) ) *K 
ZH ( I ) =SPRAD*COS ( angle ( I ) ) *K 
write ( 6 , * ) angle ( I ) ,RH(I) ,ZH(I) 

888  continue 

OPEN ( 2 , FI LE=PHI PTS ) 

OPEN ( 3 , FILE=RHOPTS ) 

OPEN ( 4 , FILE= ' etheta . m ' ) 

OPEN ( 8 , FILE= ' erho . m ' ) 

OPEN ( 9 , FILE= ' ephi . m ' ) 

OPEN(10,FILE='ang.m' ) 

READ ( 2 , * ) NPHI 

READ ( 3 , * )  NRHO 

DO  20  M=1 , NPHI 

READ ( 2 , * , END=2  0 ) XP (M) ,AP(M) 

20  CONTINUE 

DO  30  N=1 , NRHO 

READ ( 3 , * , END=30) XR (N) , AR (N) 

30  CONTINUE 

C  OUTER  LOOP  :  SOLVES  FOR  THE  SPHERICAL  COMPONENTS  OF  E  AT 
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C  EACH  POINT  ON  THE  SURFACE  OF  THE  OGIVE. 

C  MIDDLE  LOOP:  INTEGRATE  OVER  PHI . . . -PI<PHI<PI . 

C  INNER  LOOP:  INTEGRATE  OVER  RHO. . .  0  <RHO<ANT 
DO  40  1=1, NP 

Z=ZH(I)/K 
RZ=RH (I)/K 

S=RZ/SQRT (RZ**2+Z**2) 

C=  Z/SQRT (RZ**2+Z**2 ) 

DEG=angle (I) *180/PI 
sumx= ( 0 . 0 , 0 . 0 ) 
sumy= ( 0 . 0 , 0 . 0 ) 
suxnz=(0.0,0.0) 

DO  50  M=1 , NPHI 

PHIPRIME=PI*XP(M) 

DO  60  N=1 , NRHO 
RP= (ANT*XR (N) ) +ANT 

R=SQRT( (RZ-RP) **2+ (Z**2 ) + (4 *RZ*RP*SIN ( (PHI-PHIPRIME) /2 ) **2 ) ) 
Gl= ( ( (K*R) **2 ) -1- ( JK*R) )/R**3 
G2= ( 3+ (3*JK*R) - (K*R) **2)/R**5 
X1=JK* ( (RP*COS (PHIPRIME) *SIN (SCAN) ) -R) 

Yl= (RZ*COS (PHI) ) -(RP*COS (PHIPRIME) ) 

X2=Y1**2 

Y2= (RZ*SIN (PHI) )-(RP*SIN (PHIPRIME) ) 

CC=CON*ANT 

sumx=surax+ (CC*AP (M) *AR (N) * (G1+ (X2  *G2 ) ) *CEXP ( XI ) *RP) 
sumy=sumy+ ( CC*AP (M) *AR(N) *Y1*Y2*G2*CEXP (XI) *RP) 
sumz-sumz+ (CC*AP (M) *AR(N) *Y1*Z*G2*CEXP(X1) *RP) 

60  CONTINUE 
50  CONTINUE 

C  CALCULATE  THE  SPHERICAL  COMPONENTS  AT  THE  SURFACE. 

ETHETA ( I ) = ( ( C*COS ( PHI ) *sumx) + (C*SIN (PHI) *sumy) - (S*sumz ) ) 

ERHO(I)  =( (S*COS(PHI)*sumx)+(S*SIN(PHI) *sumy ) + (C*sumz ) ) 
EPHI (I)  =( (-SIN (PHI) *sumx) + ( COS ( PHI ) *sumy) ) 

ET=CABS ( ETHETA (I) ) 

ER=CABS ( ERHO ( I ) ) 

EP=CABS ( EPHI ( I ) ) 

WRITE ( 4 , * ) ET 
WRITE (8, *)ER 
WRITE (9 , * ) EP 
WRITE (10,*) DEG 
40  CONTINUE 
STOP 
END 


4 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OGIVE 

4  SEPTEMBER  1991 
R.M.  FRANCIS 
13  JANUARY  1992 


SUBROUTINE 
DATE 

PROGRAMMER 
REVISED 
COMMENTS 

THIS  SUBROUTINE  WILL  GENERATE  DATA  FOR  A  BODY  OF  REVOLUTION 
(BOR)  IN  THE  FORM  OF  AN  OGIVE. 

DIMENSIONS  ARE  NORMALIZED  TO  WAVELENGTH. 

ZH  =  Z  CO-ORDINATE  *  2*PI 
RH  =  RADIUS  *2*PI 


C* **************** 

SUBROUTINE  OGIVE (NP, ZH, RH) 

C  NP  =  NUMBER  OF  POINTS  ON  THE  OGIVE  SURFACE,  MAXIMUM  =  1000 
C  ZP  =  Z PRIME ,  THE  POSITION  ON  Z  WHERE  THE  RADIUS  OF  CURVATURE 
C  STARTS 

C  B  =  BASE  RADIUS 

C  RS  =  RADIUS  OF  CURVATURE  IN  THE  RZ , Z  PLANE  OF  THE  OGIVE 
INTEGER  I , NP 

REAL  RH (400) ,ZH(400) , ZP, BASE , RS , ZCOORD, RADIUS , PI 
PI=  3.141592654 

C  INPUT  THE  VARIABLES  FOR  THE  OGIVE, ZP, B,RS,NP 

WRITE(6,*) 'ENTER  AN  EVEN  NUMBER  OF  POINTS,  MAX  IS  400' 
READ ( 5 , * ) NP 

WRITE ( 6 , * )  'ENTER  SURFACE  CURVATURE  (wavelengths)' 

READ ( 5 , * )  RS 

WRITE  (6,*)  'ENTER  ZPRIME, WHERE  CURVATURE  STARTS  (wavelengths)  ' 
READ ( 5 , * )  ZP 

WRITE (6,*)  'ENTER  BASE  RADIUS  (wavelengths)' 

READ (5,*)  BASE 
C  PERFORM  CALCULATIONS 

ZMAX=  SQRT( (2*BASE*RS) -BASE**2)  +  ZP 
DZ=  ZMAX/ (float (NP)-l.) 

DO  10  1=1, NP 

ZCOORD=  ( float (I) -1 . ) *DZ 
ZH ( I ) =2 . * PI* ZCOORD 

IF  (ZCOORD. LE.ZP)  THEN 
RADIUS=BASE 

ELSE 

RADIUS=SQRT (RS**2- ( ZCOORD-ZP) **2 ) + ( BASE-RS ) 
ENDIF 

RH ( I ) =2 . *PI*RADIUS 
10  CONTINUE 
RETURN 
END 
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2.  Source  Code  for  TEST.M 


%  PROGRAM  :  TEST.M 

%  DATE  :  3  OCTOBER  1992 

%  REVISED  :  15  FEBRUARY  1992 

%Program  TEST.M —  generate  analytic  soln.  for  circular 
%  aperture. 

%saves  the  data  to  "mat  file",  ctest. 
npts=input (' enter  the  number  of  points:') 

R=input( 'enter  the  distance  (R):') 

r=input( 'enter  the  radius  (in  wavelengths):  ') 

scan=input( 'enter  the  scan  angle (degrees) ' ) 

sa=scan*pi/180 ; 

eta=  120*pi; 

k=2*pi ; 

j=sqrt (-1) ; 

nmax=2* (npts-l)  ; 

cl=j  *k*eta*rA2*exp (-j *k*R) /R; 

%  calculate  the  solution  from  0  to  90  degrees. ******* 
for  i=0: npts-l 
ang(i+l)=pi*i/nmax; 
deg ( i+l) =90*i/ (npts-l) ; 
fac=0 . 5*cos (ang (i+l) ) ; 
if  sin (ang (i+l) ) -sin (sa) ==  0.0 
a ( i+l) =eps; 

else 

a(i+l)=2*pi*r* (sin (ang (i+l) ) -sin(sa) ) ; 

end 


end 


z ( i+l) =fac*abs (cl*bessel ( 1 , a ( i+l) ) /a ( i+l ) ) ; 


np=npts ; 

save  ctest  deg  z  np 


3.  Data  file 

0 .  OOOOOOOe-^OO 
4 . 7368421e+00 
9 . 4736842e+00 
1 . 4210526e+01 
1 . 8947368e+01 
2. 3684211e+01 
2. 8421053 e+ 01 
3 . 3157895e+01 
3 . 78947  37e+01 
4 . 2631579e+01 
4 . 7368421e+01 
5 . 2105263e+01 
5. 6842105e+01 
6. 1578947e+01 
6 . 63 15789e+01 
7 . 1052632e+01 
7. 5789474e+01 


CTEST. MAT 

5 . 9217626e+01 
5 . 7051528e+01 
5 . 0941450e+01 
4 . 1944793e+01 
3 . 1506384e+01 
2. 1103366e+01 
1. 1933350e+01 
4 . 7301405e+00 
2 . 67 30656e-01 
3 . 2224083e+00 
4 . 5450714e+00 
4 . 7331488e+00 
4 . 2520406e+00 
3 . 4672834e+00 
2 . 6244918e+00 
1 . 8607876e+00 
1 . 2307921e+00 


20 
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8 . 05263 16e+01  7 . 3451071e-01 
8. 5263158e+01  3 . 4018017e-01 
9 . 0000000e+01  2.4512433e-16 

4.  Source  code  for  SCANPLT.M 

% PROGRAM  :  scanplt.m 

%DATE  :  15  February  1992 

%Comments  :  This  program  is  a  plotting  routine  for  the  data 
%  points  generated  by  the  test  program  CIRCSUB.F  and  the 
%  program  TEST.M.  Enter  parameters  iaw  the  runs  of 
%  TEST.M  and  CIRCSUB.F. 
load  ang.m 
load  etheta.m 
load  ctest 

subplot (211) , plot (deg, z,ang,etheta, ' . ' ) 

title ('TEST  OF  SUBROUTINE  CIRCTHETA '), xlabel (' degrees ' ) 
y label ( 1 magnitude ' ) 

text (. 1, . 4 analytical  (solid),  CIRCTHETA  (.)',' sc ' ) 

text ( . 1 , . 3 , 1  PARAMETERS : ' , ' sc ' ) 

text (. 1, .25, 'number  of  points  =  90', 'sc') 

text ( .5, .25, 'scan  angle  =  60  deg.', 'sc') 

text ( . 1 , . 2 , ' CNRHO  *  20',' SC ' ) 

text (.5, .2, 'CNPHI  =  60', 'sc') 

text (. 1, . 15, 'distance  =  2500','sc') 

text (. 5, . 15, 'aperture  radius  =  5', 'sc') 
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APPENDIX  E.  MATH  NOTES 
A.  PLOTS  GENERATED  BY  SUBROUTINE  OGIVE 

Plots  generated  by  the  subroutine  OGIVE  do  not  have  equal 
spacing  in  the  angle  0.  As  the  angle  increases  the  points  are 
spaced  more  closely.  Equation  2.1  defines  r(z).  It  follows 
that 

0  =  tan'1  |  1  j  (E.  l) 

as  z  approaches  zero,  0  approaches  90  degrees.  The  rate  at 
which  0  approaches  90  degrees  decreases  as  z  decreases  since 

- ?1 -  (E.2) 

d  z  z2  +  r(z)2 

clearly  demonstrates  in  the  limit 

limit  -f-|  =  0  (E .  3 ) 

z-  0 

Where  the  subroutine  OGIVE  generates  the  points  RH  and  ZH,  the 
density  of  plotted  values  increases  with  angle  0. 
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B.  PLOTS  GENERATED  BY  SUBROUTINE  TESTSPHERE 

In  the  subroutine  TESTSPHERE,  the  variables  RH  and  ZH  are 
defined  in  terms  of  a  sphere  radius  (SPRAD)  and  the  angle  0  . 
The  following  pseudocode  applies 

e  (i)  =  escart  ♦  a  0 

RH(i)  =  SPRAD  *  sin  ( 0  (  i  )  )  (E-4> 

ZH(i )  =  SPRAD  *  cos  (  0  (  i  )  ) 

where  A0  is  a  constant  given  by 

and  NP  is  the  number  of  points  on  the  BOR. 


A  0  =  — (E.5) 
NP  -  1 

C.  NEAR  FIELD  EQUATIONS  FOR  A  CIRCULAR  APERTURE 

The  Equations  3.13  to  3.22  were  adapted  from  Equations  (6- 
108a)  to  (6-108c)  on  page  283  of  reference  5.  To  implement  the 
near  field  integrals  for  a  circular  aperture  bounded  by  a  BOR 
the  cartesian  representation  of  the  equations  in  reference  ?? 
were  converted  to  polar  coordinates  as  follows 

X1  =  p'  COS  (J)' 

y'  =  p'  sin  <t>'  (E.  6) 

ds /  =  dx'  dy'  =  p'  dp'  d4>; 
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Where  the  primed  notation  denotes  a  source  coordinate. 
Straight  forward  substitution  of  Equation  (E.6)  into  the 
equations  of  reference  5  yeild  Equations  3.13  to  3.22. 

The  cartesian  components  of  the  electric  field  E(x,y,z) 
are  transformed  to  spherical  components  by 

Ee  [cos  (0)  cos  (<j>)  cos  (6)  sin  (<J>)  -sin  (0)1  Ex 

Ep  =  sin (0) cos  (4>)  sin (0)  sin (<{>)  cos(0)  Ey  (E*7> 

£■.  -sin(4>)  cos($)  0  e, 

▼J  2. 


Figure  E.l.  Polar  Geometry  for  Conversion  of  E(x,y,z)  to 
Spherical  Components, E(R,0,^) . 

Figure  E.l.  illustrates  the  geometry. 
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